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Abstract 

A  numerical  methodology  which  determines  the  quality  (or  robustness)  of  a- 
posteriori  error  estimators  is  described.  The  methodology  accounts  precisely  for 
the  factors  which  affect  the  quality  of  error  estimators  for  finite-element  solutions 
of  linear  elliptic  problems,  namely,  the  local  geometry  of  the  grid  and  the  structure 
of  the  solution.  The  methodology  can  be  employed  to  check  the  robustness  of  any 
estimator  for  the  complex  grids  which  are  used  in  engineering  computations. 
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1  Introduction 


A-posteriori  error  estimation  has  become  a  key  feature  of  practical  finite-element 
analysis.  Because  of  their  practical  importance  error  estimators  have  been  the 
focus  of  intensive  research;  see  for  example  [1-53]  and  the  references  in  these  pa¬ 
pers.  While  some  a-posteriori  error  estimators  have  been  analyzed  mathematically 
(e.g.  [1],  [6],  [9],  [13],  [14],  [15],  [16],  [17],  [23],  [38],  [39])  many  estimators  have 
been  derived  by  purely  heuristic  reasoning.  Usually  the  estimators  are  validated 
numerically  on  a  set  of  benchmarks  (example  problems)  which  are  selected  in  an  ad- 
hoc  manner.  Most  benchmark  computations  fail  to  isolate  the  basic  factors  which 
influence  the  performance  of  estimators  and  can  motivate  wrong  conclusions.  In 
this  paper  we  present  a  clearly  formulated  objective  validation  principle  for  error 
estimators  which  takes  into  account  the  major  factors  influencing  the  performance 
of  estimators  in  the  case  when  the  element  is  not  at  the  boundary  and  the  ex¬ 
act  solution  is  smooth  (in  the  neighborhood  of  the  element).  The  methodology  is 
completely  numerical  and  can  be  used  even  when  the  definition  of  the  estimator  is 
unknown  and  is  given  only  as  a  black-box  computer  program. 

In  practice,  we  are  interested  to  have  an  accurate  estimate  of  the  error  in  a  cell 

of  the  mesh  T/,  (we  will  use  the  term  cell  to  refer  to  a  small  patch  of  elements; 
the  cell  may  consist  of  a  few  (possibly  one)  elements).  The  performance  of  any 
error  estimator  in  u/£  depends  on  the  local  geometry  of  the  mesh  in  a  slightly  bigger 
patch  a/J  which  includes  u>£  in  its  interior  (see  Fig.  1)  and  on  the  local  structure  of 
the  solution  and  no  heuristic  benchmarks  can  properly  account  for  these  factors. 
The  methodology  given  below  enables  us  to  focus  in  the  cell  of  interest  and  to  study 
the  robustness  of  any  error  estimator  (even  if  it  is  only  available  as  a  black-box 
subroutine)  for  the  actual  geometries  of  the  grids  which  are  used  in  the  engineering 
computations.  The  methodology  requires  the  solution  of  relatively  small  problems 
in  the  regions  of  interest  and  is  inexpensive.  In  contrast,  benchmarks  require  global 
computation  (which  can  be  expensive)  and  do  not  lead  to  reliable  conclusions. 

The  quality  of  an  error  estimator  in  the  cell  u/q  is  measured  by  the  effectivity 
index 


»  S4:=={  £ 
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where  |||ej,|||ll,k  is  the  norm  (of  interest)  of  the  error  over  u>£,  is  an  error  esti¬ 
mator  for  this  norm  which  is  computed  in  terms  of  element  error-indicators  rjr,  r 
denotes  an  element  in  the  mesh  7*.  In  this  paper,  we  will  consider  estimators  for 
the  energy-norm  of  the  error.  The  methodology  of  the  paper  can  also  be  used  to 
study  the  quality  of  estimators  for  other  norms.  Let  ft  denote  the  domain  of  the 
problem;  in  [1],  [6],  [9],  [13],  [14]  it  has  been  shown  that  the  range  of  the  global 
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effectivity  index,  «n,  exists  for  several  estimators  based  on  residuals,  namely  there 
exist  constants  0  <  C2  <  C§  <  +oo  such  that 

0<Cg<Kn<Cg<oo  (2a) 


The  two-sided  inequality  (2a)  has  been  proven  under  very  general  assumptions 
about  the  exact  solution  (it  is  only  required  that  the  exact  solution  has  finite- 
energy),  reasonable  assumptions  about  the  regularity  of  the  data  (all  practical 
cases  are  included)  and  under  mild  restrictions  on  the  regularity  of  the  mesh  (see 
the  details  in  [6],  [9],  [13]).  Inequality  (2a)  can  also  be  written  in  the  form 

^I&<ll|e‘lllo<gff6.  (26) 

which  expresses  an  equivalence  between  the  global  norm  of  the  error  and  the  estima¬ 
tor.  Practical  values  for  the  equivalency  constants  C2  and  cannot  be  obtained 
for  a  given  estimator  unless  further  information  is  known  about  the  class  of  solu¬ 
tions  of  interest  and  the  finite-element  meshes  employed.  A  concrete  example  of 
how  the  constants  can  be  estimated  in  the  case  of  a  simple  residual  estimator  was 
given  in  [50]  and  [51].  The  values  of  and  Cf)  depend  strongly  on  the  geometry 
of  the  grid  and  (relatively  weakly)  on  the  smoothness  of  the  solution;  the  geometry 
of  the  grid  must  be  understood  in  connection  with  the  differential  operator  (see 
[50]  for  the  details). 

It  can  also  be  shown  (see  [53]  and  the  outline  given  below)  that  under  reasonable 
assumptions  about  the  grid  we  can  determine  the  asymptotic  range  of  the  effectivity 
index  for  any  estimator  in  any  small  interior-cell  u (a  cell  which  is  separated 
from  the  boundary  by  several  mesh-layers)  of  the  grid;  i.e.  there  exist  constants 
0  <  C£°  <  Cp®  <  +oo  which  depend  only  on  the  local  geometry  of  the  grid  in  u/J 
and  a  few  mesh-layers  around  it  (the  geometry  of  the  grid  in  a  sufficiently  large 
patch  u;*  which  includes  in  its  interior)  such  that  (as  the  local  mesh-size  in 
and  w*  tends  to  zero) 

0  <  cf  <  kuh  <  Cf  <  oo  (3) 

The  constants  C£° ,  C'y  are  the  best  possible  constants  (i.e.  they  can  be 
achieved)  over  an  entire  class  of  smooth  solutions  which  occur  in  the  field  of  appli¬ 
cation.  When  the  values  of  the  constants  C£® ,  C£®  are  known  we  can  also  manu¬ 
facture  an  upper-( resp.  lower-)  estimator  version  of  £uh  denoted  by  (resp.  £"k) 
defined  by 


(4a) 
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We  then  have 


(46) 


The  lower-estimator  version  of  an  estimator  may  be  employed  to  drive  adaptive- 
refinement  algorithms  in  order  to  avoid  over-refinement  while  the  upper-estimator 
version  may  be  used  to  guarantee  a  safe  stopping-criterion. 

Let  us  define  the  robustness  index  Huh  (0  <  7luk  <  oo)  which  expresses  the 
reliability  of  the  estimator: 

*„}  :=  max((|l  -  Ctf  I  +  U  -  I)  ,  (U  -  il  +  H  ~  *4 I)}  (5«) 

Cv  '-'L 

The  robustness  index  expresses  the  deviation  of  k  and  —  (see  (2a),  (2b))  from  the 

ideal  value  n  =  1.  (Hence  =  0  is  the  ideal  value  for  the  robustness  index.) 
The  robustness  of  an  estimator  for  a  given  class  of  meshes  T  =  {Th}  is  given  by 

tt(r)  :=  (5fe) 

Here  C(Tk)  is  the  set  of  interior  cells  in  the  mesh  7/,.  The  robustness  index  defined 
in  (5b)  characterizes  the  performance  of  the  estimator  for  elements  in  the  interior 
of  the  domain.  Because  most  elements  are  in  the  interior  of  the  domain  (where 
the  solution  is  smooth)  the  validation  of  the  estimators  has  to  be  made  especially 
for  the  case  presented  in  this  paper.  The  validation  of  the  estimators  for  elements 
located  at  the  boundary  and  for  unsmooth  solutions  will  be  given  in  subsequent 
papers. 

If  an  estimator  does  not  display  reasonable  robustness  for  the  interior-cells, 
i.e.  if  72.(7”)  is  too  large,  the  estimator  is  unreliable  and  should  not  be  considered 
any  further.  The  robustness  index  depends  on  7”,  the  family  of  meshes  under 
consideration.  (Hence  restrictions  placed  on  the  mesh-generator  could  possibly 
increase  the  reliability  of  an  estimator.) 

The  robustness  index  is  an  objective  quantitative  characterization  of  the  relia¬ 
bility  of  an  estimator.  Hence,  analogously  as  the  effectivity  index,  the  robustness 
index  of  an  estimator  should  be  reported. 
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Following  this  Introduction  we  present  the  definition  of  the  model  problems 
(linear  elasticity  and  heat-conduction),  we  describe  two  types  of  error  estimators, 
we  present  the  methodology  for  the  computation  of  the  robustness  index  and  we 
outline  its  theoretical  justification.  Finally,  we  give  examples  which  demonstrate 
how  the  robustness  index  of  error  estimators  for  complex  finite-element  meshes  can 
be  computed  and  how  it  is  possible  to  increase  the  reliability  of  an  estimator  by 
proper  selection  of  its  various  parameters. 


2  The  model  problems 

We  shall  consider  the  vector-valued  boundary-value  problem 

Ii(u)  :=  -  XI  Di  (<MU))  =  fi  “  ft 
j= i 


tij  =  0  on  I'd 

2 

X3  =  &  on  IV 

j= 1 


(6) 


where  i  =  1,  2. 

Here  12  C  R2  is  a  bounded  domain  with  boundary  dCl  =  Ip  U  IV ; 

n  :=  (ni,n2)  is  the  outward  pointing  unit-normal  on  IV; 

fi,  i  =  1,2  are  the  components  of  the  load- vector  ( body- force ); 

U,  i  —  1,2  are  the  components  of  the  normal-flux  vector  (traction)  applied  on  IV; 
IV  ^  0,  rDnrw  =  0;u=  («i,u2)  is  the  solution-vector  ( displacement ); 

6jj(u)  ;=  -{DjVi  +  DiVj)  ,  i,j  =  1,2  (7a) 

2 

<7*j(u)  :=  £  Oijueui u)  ,  i,j  =1,2  (76) 

M=i 

are  the  components  of  the  flux  ( strain ,  stress); 

Oijki ,  t,j,k,l  =  1,2,  are  the  material-coefficients  ( elastic  constants)  which  satisfy 

Uijki  —  <*jUk  =  aikij ,  i,j,k,l=  1,2  (8a) 
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3  2 

E  “vfW  >  c  52  «<*  »  c  >  0 ,  V  ey  =  e*  (86) 

*j'=i 


(Conditions  (8a),  (8b)  are  satisfied  for  linear  elasticity;  in  the  case  of  isotropic 
plane  elasticity  we  have  Oijiu  =  +  Su&kj)  +  A 8*6 jt  where  is  Kronecker’s 

delta  and  X,p  are  Lame’s  constants.) 

We  also  introduce  the  scalar  elliptic  boundary-value  problem  (heat-conduction 
in  orthotropic  medium ),  namely 


L\u)  :=-  E  Dk(KklDlu)  =  f 

in  H 

M=i 

u  =  0 

on  I'd 

2 

E  9*(u)n*  =  * 

fc=l 

on  VN 

(9) 


Here  qh( u)  :=  £  KuDtu ,  6  =  1,2  are  the  components  of  the  flux-vector  (heat- 

t=i 

flux)  and  Ku,  6,/  =  1,2,  are  the  entries  of  the  thermal-conductivity  matrix  which 
satisfy 


Ku  =  Ka,  ,  k,l  =  1,2 


(10a) 


0  <  AWtf  +  $)  <  E  <  *»«(£  +  €) ,  £  €  R*  (106) 

fc/sl 

Here  denote  the  principal  values  of  the  thermal-conductivity  matrix. 

Let  us  now  cast  the  model  problem  in  variational  form.  Let  us  denote 

H1  :=  |  (t>i,vj) :  v,  €  #*(0)  J 
Hj>0  :=  |  (t>x,V2) :  Vi  €  #*(8) ,  t><  =  0  on  TD  J 

iMii4»:=  (emu)  » ivii^:=  (emu) 
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with  ||v«||i,fl  (|i>j|i,n)  being  the  usual  Sobolev  norm  (seminorm).  The  vari¬ 

ational  form  of  the  boundary-value  problem  (6)  is  now  posed  as: 

Find  u  €  such  that 

Bn(u,v)  =  f  £>.  +  /  v  v  €  H£d  (11) 

Ja  i= i  Jrtf  i=  i 

where  the  bilinear  form  Bn  :  Hrc  x  HpD  — +  R  is  defined  by 

t 

Bn(u,v)  :=  /  D  “HuDtUjDkVi .  (12) 

•,n 

The  energy-norm  over  any  subdomain  5  C  12  is  defined  by 

HMDs  :=  y/Bs(v,v)  (13) 

where  Bs( u,v)  has  the  obvious  meaning. 

In  the  case  of  the  scalar  elliptic  problem  given  by  (9)  the  bilinear  form  is  given 

by  B'a(u,v)  :=  /  D  KhiDt\iDhv .  The  weak-solntion  of  (9)  satisfies: 

Jn  m=i 

Find  u  6  Hj-d  :=  jv  6  /^(fl)  :  v  =  0  on  I'd  j  such  that 

B'a(u,v)  =  j[/»  +  V  *  6  fl*,  (14) 

The  energy-norm  in  any  subdomain  S  C  12  is  defined  by  |||v|||5  :=  ^B's(v,v)  . 

Let  T  =  {Th}  be  a  family  of  meshes  of  triangles  or  quadrilaterals  with  straight 
edges.  It  is  assumed  that  the  family  is  regular  (i.e.:  for  the  triangles  the  minimal 
angle  of  all  the  triangles  is  bounded  below  by  a  positive  constant,  the  same  for 
all  the  meshes;  for  the  quadrilaterals  see  conditions  (37.40)  in  Ciarlet  [58]).  The 
meshes  are  not  assumed  to  be  quasiuniform.  Let  us  introduce  the  finite-element 
spaces 


Sph  :=  {u  €  H1  :  €  Vp[t),  i  =  1,2,  fc  =  l,...,M(Th)}  (15) 

where  is  the  mapping  function  for  the  &th  finite-element  which  maps  either  a 
standard  triangular  element  (using  an  affine  transformation)  or  a  standard  quadri¬ 
lateral  element  (using  a  bilinear  transformation)  onto  the  kih  finite  dement,  f 
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denotes  a  standard  element,  M(Th)  is  the  number  of  elements  in  the  mesh  7J,, 
V,(f)  denotes  the  set  o'  polynomials  of  degree  p  over  ?.  We  let 

Shx„  :=  n  . 

The  finite  element  solution  u*  (for  the  elasticity  problem)  satisfies: 

Find  u*  €  S£  rp  such  that 

Bo(«‘,v*)  = /  f /,<,?  +  /  f*”*  V  v*  €  s;  r  (16) 

The  error  is  e*  :=  u  —  u* . 

Remark  2.1.  We  addressed  only  (for  simplicity)  the  model  problems  for  which 
the  differential  operator  L  (or  L')  is  homogeneous  with  constant  coefficients.  The 
theory  and  the  procedure  holds  for  the  general  case,  for  non-homogeneous  operators 
with  non-constant  (but  smooth)  coefficients. 

In  the  following  we  give  two  representative  classes  of  estimators  for  the  energy- 
norm  of  the  error  and  we  employ  the  model  methodology  to  check  the  robustness 
of  various  versions  of  estimators  from  these  clashes.  Below  we  define  the  estimators 
for  the  elasticity  problem  (the  estimators  for  the  scalar  model  problem  (9)  may  be 
obtained  by  analogy;  see  also  [53]). 

3  Element-residual  error  estimators 

3.1  Implicit  element-residual  estimator 

We  introduce  notations  needed  for  the  description  of  the  estimators  (see  also  [8], 
[14]).  For  each  closed  triangle  (quadrilateral)  r  €  7*,  we  denote  by  E(r)  and  N(r) 
the  set  of  its  edges  and  vertices,  respectively.  We  define  the  local  bubble-spaces  of 
polynomials  (see  [8],  [11],  [12],  [14],  [21],  [22],  [23]) 

H,  :=  jw  :  Wi  €  Vp+i ,  «\(Jf)  =0  V  X  €  JV(r) ,  *  =  1,2  j  , 

over  each  element  r  €  7*. 

We  define  the  interior-residual  in  element  r  as 

rT  :=  — L(u*)+f  (17) 

where  L(n*)  :=  (Li(u*),La(u*))  ,  f:=  (/i,/a), 
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and  the  jump  or  edge-residual  associated  with  the  edge  e 

J«:=  [ff(ufc|TMt)-ff(ufc|^,)]n  (18) 

where  n  is  the  unit-normal  for  the  edge  e  and  and  are  defined  as  in  Fig.  2. 
The  residual-functional  for  element  r  is 


«7>(v):=  /vtt  +  ^  £  v€Hx(r) 


(19) 


<e*W 


We  now  define  the  element-residual  estimator  for  the  model  problem  (6)  (we 
give  the  estimators  only  for  elements  in  the  interior  of  the  domain). 

The  element  error  indicators  for  the  implicit  element-residual  estimators  for 
elements  of  any  degree  p  are: 


Vr  :=  Hkrlllr  (20) 

where 


er  €  Hr  :  BT(er,vr)  =  /y(vT)  ,  V  vT  €  Hr  (21) 


3.2  Equilibrium  of  the  residuals 

The  residual  data  for  the  local  problems  (21)  are  equilibrated  if  the  following 
consistency  conditions  ( equilibrium  equations)  are  satisfied, 


jf T*  +  \  £  /J‘ =  0il  +  0ij  | 


«€J?(r) 


/*xr r  +  J  E  /xx  J,=  0i8 


«€*(r) 


(22) 


Here  ii,  ij  denote  the  unit- vectors  along  the  global  coordinate  directions  in  R3  and 
is  the  out-of-plane  unit-vector.  The  element  equilibrium  equations  (22)  will  hold 


if 


E{^r(«iJVi)ii+^,(ij«i)ii}  =0i,  +0i, 
1=1 

nv 

x  {^.(iiAT^ii  +  ^(isWOh}  =0i, 

»=i 


(23) 
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It  is  assumed  that  the  set  of  functions  {#»}£?!  satisfies  /V*  =  1 ,  V  r  6  7J, . 

t=i 

The  equilibrium  conditions  for  element  r  will  be  satisfied  automatically  if 

•7>(i iNi)  =  0,  £  =  1,2,  *=l,...,nv  (24) 

These  conditions,  however,  are  not  expected  to  hold  for  general  meshes  and  solu¬ 
tions;  below  we  discuss  ways  of  correcting  the  definition  of  Tr  to  satisfy  conditions 
(24). 

Remark  3.1.  A  possible  choice  for  the  functions  {N}^  is  a  set  of  Lagrangian 
basis  functions  of  degree  q  <  p  (p  is  the  polynomial  degree  of  the  finite-element 
space)  defined  over  7J,.  In  the  discussion  below  we  let  9=1,  Nil  is  the  linear  (for 
triangles)  or  bilinear  (for  quadrilaterals)  element-shape-function  which  corresponds 
to  the  t-th  vertex  and  nv  is  the  number  of  vertices  for  element  r  (nv  =  3  for 
triangles,  nv  =  4  for  quadrilaterals). 

3.2.1  Ladeveze’s  flux-splitting  technique 

The  definition  of  in  (19)  may  to  be  modified  in  order  to  satisfy  (24)  by  letting 

*T*(v)  :=  Triy)  +  /  V  •  eT  ,  v  €  Hl(r)  (25) 

J&r 

Here,  0T  €  (lJ(3r))2  is  the  correction  of  the  edge-residuals  for  the  element  r. 

We  let 


«M.  =  C(e>;+eta)  ,  .6*( r)  (26) 

where, 


f+1,  = 

C  = 

v  -1,  if  r  =  r^, 

where  it  is  assumed  that  the  edge  normal  n  has  been  assigned  to  the  edge  c  in  an 
arbitrary  but  unique  way  and  rta  and  tm  are  defined  as  shown  in  Fig.  2.  Here,  the 
linear  functions  ipl  are  defined  as, 

«  -  A  (2A;  -  a;+1)  ,  *=1,2  (27) 

where  A*,  k  =  1,2  are  the  linear  shape-functions  defined  over  the  edge  c  and  the 
subscript  k  + 1  is  taken  modulo-2.  Using  this  definition  of  0T|,  we  can  decouple  the 
problem  of  determination  of  0T|t  for  the  whole  domain  into  small  local  problems 
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involving  only  patches  of  elements  connected  to  node  X ,  as  shown  in  Fig.  3a.  We 
have 


ef:=^0T  |,AJf  k  =  1,2,  eeE{r)  (28) 

Thus,  for  each  patch  around  a  node  X  we  obtain  a  linear  system  (see  [24],  [8],  [32], 
[59]  and  [60]  for  the  details) 

!^x( *  ®r*)  <i>x  =  r(ll<f>X )  ,  /  =  1, 2  fc  =  1, . . . , Nx  (29) 

where  <f>x  denotes  the  elementwise  affine  basis  function  (for  meshes  of  triangles) 
which  corresponds  to  the  vertex  X,  rx  denotes  the  k- th  element  connected  to  the 
vertex  X,  Nx  is  the  number  of  elements  (or  edges)  connected  to  the  vertex  X. 

The  procedure  outlined  above  has  been  developed  by  Ladeveze  [59].  The  linear 
system  (29)  has  a  one  parameter  family  of  solutions.  Specific  choices  of  solutions 
are  suggested  in  [8],  [25],  [32],  [60]  and  [61];  in  the  numerical  implementations 
we  employed  these  choices.  Below  we  give  the  definition  of  the  edgewise-linear 
correction  8r  (g-degree  (q  >  1)  polynomial  representation  of  the  correction  8T  over 
each  edge  is  also  possible)  which  results  by  using  the  equilibration  procedures  in  [8], 

[25] ,  [32],  [60]  and  [61].  Let  us  consider  the  interior- vertex  X  and  let  us  also  denote 
ti,  i  =  1, . . . ,  Nx  the  edges  connected  to  X.  We  will  determine  the  coefficients 

which  is  associated  with  the  edge  e  and  the  vertex  X  and  is  employed  in 

(26) .  Here  the  index  v(e,X)  identifies  the  local  enumeration  of  node  X,  as  used  in 
(26)  for  the  unknowns  associated  with  the  edge;  see  Fig.  3b. 

a.  Bank  and  Weiser’s  equilibration: 

A  solution  of  (29)  can  be  obtained  by  letting 

<=1,2 

i, .  e?"-*'  =  i,  ■  e£;- Jr>  -  r^u+x) , 

This  choice  of  solution  has  been  suggested  in  [8]. 

b.  Ladeveze  *s  equilibration : 

Here  the  coefficients  8^*^  =  £  (L  •  8rf*”^)i*,  i  =  1, . . . ,  N  are  selected  such 

i=i 

that 


(30) 


l  —  1,2;  »  =  2,...,JV 


«=i 


(31) 
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is  minimized  for  l  =  1  and  1  =  2  separately,  where  Wi  is  the  weight  associated  with 
each  edge  e*.  In  [25]  and  [32]  the  weights  to*  have  been  taken  as  1.  For  this  choice 
of  weights,  the  solution  of  (24)  is  given  explicitly  as, 


U  ■  =  77  £(*-•  +  IWfWl)  ,  l  =  1,2 

»=1 


ii  -vi''*  -rTx(n*x),  1  =  1,2 

i, -Vi**'  =i,-«2«-*) ,  <  =  1,2;  i  =  2,3, 


>  (32) 


If  we  take  the  weights  Wi  :=  |e^|  1 ,  as  suggested  in  [60]  and  [61],  we  can  obtain  a 
different  8r  satisfying  (29). 


Remark  3.2.  From  eq.  (29)  we  can  see  that  the  functional  Jr^Q  with  any  of  the 
above  choices  of  the  correction  8r  as  the  above  will  satisfy  (24)  and  therefore  the 
equilibrium  conditions  (22). 


3.2.2  Ohtsubo-Kitamura  flux-splitting 

In  [29-31]  a  different  procedure  for  determining  equilibrated  force-fields  for 
an  element  r  is  described.  In  this  procedure  each  element  can  be  equilibrated 
separately  from  its  neighbors  i.e.  the  equilibrations  are  done  at  the  element-level. 
The  technique  has  as  follows: 

Let  us  define  the  residual-functional 


•*v(v)  =  / rT-v+  £  / a;j€ 

Jr  f  €S(f ) 


V  v  €  Hl(r) 


(33) 


Here  a,  is  the  geometric  splitting-factor  for  the  initial  allocation  of  the  jump  J( 
defined  by 


a?”  := 


d?" 


dp*  -I-  dp—*  ’ 


aT—*  := 


dp- 


dp*  +  dp—* 


where  dp*,  dp—*  is  the  distance  of  the  centroid  of  the  element  rto,  t**,  respectively, 
from  the  common  edge  e  (the  choice  of  the  splitting-factors  is  motivated  from  the 
one-dimensional  case  where  a  unique  flux-splitting  exists;  see  [25],  [26]  and  page 
296  in  [62]).  Following  [29-31]  we  will  assume  that  the  finite-element  solution 
has  been  computed  using  bilinear  quadrilaterals  (or  linear  triangles)  and  the  local 
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element-residual  problem  is  discretized  using  8-node  serendipity  quadrilateral  (6- 
node  Lagrangian  triangle).  Then  the  discrete  nodal-forces  for  the  local  problem 
are 

fe=  £  fc&M  (34) 


where  ,  t  =  1, . .  .,nv'  (nv1  =  8  for  quadrilaterals,  nv1  =  6  for  triangles)  are  the 
Lagrangian  shape-functions  used  in  the  discretization  of  the  local  residual-problem. 

The  consistency  conditions  for  the  discrete  element-residual  problem  in  element 
r  are 


.  nv' 

I  rT  +  £&  =  Oil  +0ij 

Jr  >=i 


f 

/  xx  rr  +  5]xi  x  &  =  0is 

Jr  i=i 


(35) 


Here  g,-  t  =  l,...,nv'  are  the  consistent  nodal-forces  which  satisfy  (35).  In  gen¬ 
eral  the  loads  &  given  by  (34)  (from  the  initial  flux-splitting)  do  not  satisfy  the 
consistency  conditions  (35).  We  seek  corrective-loads  A&  such  that 


g<  =  ft  +  Ag,  (36) 

Following  [31]  we  will  select  such  that  it  satisfies  (35)  and  has  minimal 

Euclidean-norm  in  Rnr  (this  choice  is  motivated  by  the  idea  that  the  correction 
should  be  as  small  as  possible;  see  [31]). 


4  Error  estimators  based  on  smoothening 
techniques 

Error  estimators  based  on  recovery  methods  are  given  in  [33-44].  Here  we  focus  on 
the  superconvergent  patch-recovery  technique  [40-44].  The  element  error-indicators 
for  elements  of  any  degree  p  are 


ifc.:=||<r"-<r(u*)|U., 


(37) 


where 


jr  aiju.  au 


(38) 
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where  a~^  are  the  entries  of  the  compliance  tensor,  9*  is  the  recovered  flux. 

Let  u>x  :=  U  T>  denote  the  patch  of  elements  connected  to  vertex  X.  For 
X€JV(f') 

each  patch  u>x  we  recover  the  patch-projection  ax  by  solving  the  following  least- 
squares  problem: 


inf 


||9(ufc)  -  <rT|b(«x).a-‘.{y<>2.*? 


=  IK"*)  (39) 


where  {y/}^  denotes  a  set  of  sampling  points  in  ux  and 


MU 


n#p 


n.p  :=  V 
t-i  *—• 


m=l 


E  aij(ym)  <HjU  <^w(ym) 


(40) 


9*  is  obtained  by  averaging  the  patch-projections  ox  over  the  elements  (see  [40], 

[41]). 

Note  that  the  class  of  estimators  given  in  [40-41],  [43-44]  use  only  the  approx¬ 
imate  solution  u*  to  estimate  the  error.  It  is  possible  to  modify  the  definition  of 
the  estimators  to  include  the  information  given  by  the  data  and  the  differential 
equation.  In  [42]  a  modified  patch-projection  Bx  was  defined  to  take  into  account 
the  additional  information,  namely 


F(5*)=  inf  F(«*)  (41) 

wij  €*>(~x) 

ijmia 

where 


F(o) ||o(ufc)  - *IIl»(*x),o-« +  £  H/< +  i  (42) 


i=i 


;=i 


and 


ntp 

:=  E(v(ym))  »  veC0(ux)  (43) 

m=l 

In  the  majority  of  the  practical  problems  the  body-force  vanishes  identically 
(fi  =  0,  *  =  1,2).  To  take  this  information  into  account  we  define  the  space  of 
“harmonic”  polynomial  fluxes: 

S'(wjt)  :=  U  :  €  Vp{wx) ,  *,  j  =  1,2  and  £  Dfaij)  =  0 ,  »  =  1,2}  (44) 
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We  may  now  modify  the  minimization  problem  in  (39)  by  seeking  the  minimum 
over  Sp(u>x) ;  <r*  in  this  case  will  be  called  “harmonic”  patch-projection. 

Remark  4.1.  Although  the  intention  is  to  use  superconvergence  points,  the 
estimator  performs  very  well  (as  will  be  seen  later)  also  if  the  sampling  points  are 
not  superconvergence  points  (which  for  general  meshes  do  not  exist).  In  general, 
the  performance  of  the  estimator  depends  on  the  selection  of  the  sampling  points. 
The  sampling  points  in  each  element  are  obtained  by  mapping  the  points  defined 
in  the  standard-element  for  each  element-type  (see  Fig.  4;  the  performance  of  the 
estimators  for  the  various  choices  of  sampling  points  will  be  addressed  below)  and 
may  be  selected,  more  or  less,  arbitrarily  but  must  be  invariant  with  respect  to 
cyclic  permutations  of  the  vertices  of  the  elements  (the  error  estimator  must  be 
independent  from  the  enumeration  of  the  grid). 

5  Theoretical  setting  of  the  methodology 

The  theoretical  justification  of  the  numerical  methodology  for  checking  error 
estimators  is  outlined  below  (see  [53]  for  the  complete  presentation).  The  mathe¬ 
matical  framework  is  outlined  for  special  types  of  grids  which  are  locally-periodic  in 
the  regions  of  interest  (here  we  assume  triangular  grids;  the  theoretical  setting  for 
meshes  of  quadrilaterals  with  straight-edges  follows  similar  steps).  These  grids  are 
made  locally  by  the  periodic  repetition  of  a  super-patch  as  shown  in  Fig.  5b.  The 
theory  has  asymptotic  character  (it  holds  for  sufficiently  small  size  of  the  super- 
patch)  and  can  be  also  used  to  study  the  local  asymptotic  properties  of  general 
grids  like  the  ones  shown  in  Figs.  1  and  6.  These  are  typical  meshes  produced 
by  automatic  mesh  generators  and  include  various  patches  of  different  geometrical 
structure  (The  grid  of  triangles  shown  in  Fig.  6  was  constructed  using  a  com¬ 
mercial  mesh-generator;  the  grid  of  quadrilaterals  shown  in  Fig.  1  was  obtained 
by  convening  the  mesh  of  triangles  shown  in  Fig.  6  into  a  mesh  of  quadrilater¬ 
als.  Note  that  the  grids  are  refined  in  the  interior  because  they  were  originally 
constructed  to  include  an  interior  material  interface.  We  adopted  these  grids  be¬ 
cause  they  are  rich  in  patterns  and  we  considered  model  problems  in  homogeneous 
continua.).  Here,  for  simplicity,  we  will  address  only  the  case  where  u  satisfies 
Poisson’s  equation  (— Au  =  /;  Ku  =  &iu  in  (9))  and  is  approximated  by  piecewise- 
linear  finite-elements,  (p  =  1).  Nevertheless,  the  results  could  be  easily  generalized 
to  the  general  setting  of  elasticity  in  orthotropic  medium  (see  eqs.  (6)-(8)),  to 
meshes  of  quadrilaterals  and  to  higher-order  elements. 

It  will  be  shown  that  the  robustness  index  of  the  shaded  cell  shown  in  Fig.  5c 
depends  practically  (in  the  range  of  engineering  accuracy)  on  the  mesh  in  its  neigh¬ 
borhood  (the  grid-geometry  in  the  patch  included  within  the  thick-line  peri  gram 
in  Fig.  5c)  and  not  on  the  further  extension  of  the  mesh  into  the  periodic  super¬ 
patch.  Let  us  now  specify  more  precisely  the  class  of  meshes  employed  by  the 
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mathematical  analysis.  Let  0  <  Hz0  <  Ha e  <  H°,  x°  =  (x° ,  x§)  €  H, 

S(x“,  H)  :=  {*  =  (*„ n)  :  |*i  —  *®|  <  H,  i  =  1, 2}  (45) 

and  assume  that  H°  is  sufficiently  small  that  5(x°,/f°)  C  Cl.  By  5  (x°,  H)  we 
denote  the  domain  of  S  (x°,  H)  with  its  boundary.  We  have 

Z0  :=  S(x\Hz.)  C  n,  :=  S(x°,ffn,)  • 

Further,  let  (i,j)  €  7  be  a  set  of  multi-indices,  x^  =  (x[‘J\x2  G  Cl  and 

c(x<*>, h)  :=  5(xM, h)  C  5(x°,  H°)  ,  (i,i)  G  7  (46) 

be  the  set  of  the  6-super-patches  (briefly  super-patches)  which  cover  S  (x°,  H°) 
such  that 


U  c(x^\h)  =  S(x°,H°)  (47 a) 

Me-r 

c(x(-),fc)nc(x(-),/,)  =  0  for  (iuji)?(i7j2)  .  (476) 

We  will  assume  that  S  (xP,HzQ)  (and  S(x°,lfn,))  can  be  partitioned  into  the 
set  of  super-patches 

S(x°,Hz,)=  |J  . 

(An  example  of  a  general  domain  Cl  with  the  subdomains  Z0,  Clo  and  5(x°,H°) 
is  given  in  Fig.  5a) 

Denoting  by 

c:={(*i,xa):  I  *1  |<  1,  |  x2  |<  lj  (48) 

the  unit-  (master-)  super-patch  c,  the  6-cell  is  an  6-scaled  and  translated  master 
super-patch. 

Let  farther  T  be  a  triangular  mesh  on  the  master-super-patch  (the  master- 
mesh)  and  be  the  mesh  on  c{x^\  6)  which  is  the  scaled  and  translated  image 
of  T.  Let  us  now  consider  the  family  of  regular  triangular  meshes  %  =  T  (Cl,  h)  on 
ft  with  elements  r  and  their  restrictions  7*  (x°,  H°)  on  S  (x°,  J7°).  We  will  denote 
by  the  restriction  of  7*  (x°,//°)  on  c(x^*J\6).  We  assume  that  the  meshes 
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under  consideration  are  such  that  =  7$*^  and  are  translation-invariant  in 
5  (x°,  H°).  An  example  of  such  a  mesh  is  shown  in  Fig.  5b  where  the  master  mesh 
is  shown  in  Fig.  5c.  We  note  that  we  are  dealing  with  a  one-parameter  family  of 
meshes  T(Cl ,  h).  The  parameter  2 h  has  the  meaning  of  the  length  of  the  side  of  the 
super-patch  in  5  (x°,ff°).  In  0  —  S  (x°,  H®)  it  has,  in  general,  nothing  common 
with  the  size  of  the  elements  used. 

Given  a  function  u  and  the  multi-index  a  :=  (ati,a2)  with  |a|  :=  ax  +  a2  we 
denote 


Dau  := 


dx° 1  dx? 


(49a) 


(D*u)(x)  :=  [(  52  l^attl2)(x)|  *  >  k>0,  integer  (495) 

lVM=*  j  1 

We  will  make  the  following  assumptions  about  the  exact  solution  u: 
Assumption  I 
On  S  (x° ,  Ufio ) 


\Dau\  <  K  <  oo,  0  <  |a|  <  3  (50) 

Remark  5.1.  Assumption  I  states  that  the  solution  is  smooth  in  the  neighborhood 
of  the  subdomain  of  interest  (the  subdomain  must  be  sufficiently  far  from  singular 
points,  boundaries  and  material-interfaces). 

Assumption  II 


R’=  £  [(D“«)(x")]!>0  (51) 

|«|=2 

We  will  assume  that  H°  is  sufficiently  small  depending  on  K  and  R. 

Remark  5.2.  Assumption  II  implies  that  the  principal  part  of  the  error  is  related 
to  the  (nonzero)  second-derivative  of  the  exact  solution. 

Assumption  III 

On  5(x°,Rz0),  Hz<>  <  <  H° 

<  ChfiHZo  (52) 

16 

with  0  >  —  and  where  C  is  independent  of  7(0,5),  Hz„  but  it  depends  on  K 
and  R.  Here  ||  •  ||  denotes  the  usual  L -norm. 
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Remark  5.3.  We  do  not  assume  that  u  is  smooth  in  0  outside  of  S  (x°,  H^). 
For  example,  Q  can  have  a  boundary  with  corners  (as  in  Fig.  5a)  and  hence  u  can 
be  unsmooth  in  the  neighborhood  of  these  corners.  Nevertheless  assumption  III 
makes  an  implicit  requirement  on  the  (refinement  of  the)  mesh  in  the  neighborhood 
of  these  corners  (so  that  there  is  no  significant  influence  (pollution)  of  the  errors  in 
the  neighborhood  of  the  corner  on  the  errors  in  S(x°,  Hz0) ;  for  further  discussion 
about  pollution  see  [55]). 

In  the  following  we  give  the  major  theorems  which  justify  the  methodology  of 
the  paper  (for  a  complete  presentation  with  proofs  see  [53]). 

Let  Q  be  a  quadratic  Taylor-expansion  of  u  about  x°  defined  on  S  (x°,  Hn 0 ) 
and  let  Q™T  be  its  linear  interpolant  on  the  mesh  (x°,  H^).  We  define  the 
interpolation  error 

p-=Q-QT  (53) 

which  is  a  periodic  function  on  S  (x0,/fno)  with  period  2 h  (see  [53]).  Let 

Hptn (c(x^‘J\/i))  :=  ju  €  H1  ^c(x^,J^,/i)j  :  u  satisfies  (54b) |  (54a) 

where 

u(x(1,j)  +  h,x  2)  =  i*(x(1’j)  -  h,x  2) ,  |x2  -  x£j)|  <A  | 

(545) 

«(*1  »*2  j)  +  5)  =  ^(x^x^  -  h),  |xi  -xS’j)|  <  h  J 

is  the  periodicity-condition  in  c(x^\h)  and 

~  {«  €  /?;„(c(xW»,A))  :  „[,€  7>„  r  6  (55) 

We  introduce  the  periodic  finite-element  problem : 

Find  €  S£PBR  (c(x^\hfj  Buch  that 

K> vh)  =  {P, vfc)  V  vh  €  SlPBJl  (c (x(<J),  hj)  (56) 

and 

/  iP  -  *h)  =  0  (57) 

Je  (*•>), h)  ' 
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The  function  z£  v’U  be  called  the  periodic  finite  element  approximation  of  the 
interpolation-error  p.  Let  us  remark  that  by  scaling  and  translating  the  mesh  we 
can  solve  the  finite  element  problem  in  the  master  super-patch  using  the  master 
mesh.  We  will  denote  by  z£  the  periodic  extension  of  z£  on  5(x°,  Hq0).  Let  us  now 
define  the  error  in  the  periodic  finite-element  approximation  of  the  interpolation- 
error 


4  :=  p-z'h  (58) 

The  main  theorem  of  the  paper  states  that  eh  (the  error  in  the  finite-element 
approximation)  has  the  same  asymptotic  behavior  as  rff  (the  error  in  the  periodic 
finite-element  approximation  of  the  interpolation-error  in  the  quadratic  Taylor- 
expansion  of  the  solution  about  x°)  in  the  neighborhood  of  x°. 

Theorem  1.  Let  Hz0  <  Hq0  <  H°,  and  let  assumptions  I-III  hold  and 

C,Bl  <h<  C,B%,  a  =  |  .  (59) 


Thai  if  p  = 


in  assumption  III 


|||e'*|||s(*°.ir«0)  =  ||hHI|s(x«.V*)(l  +  (60) 


where  C  is  independent  of  h  and  Hq.  ,  7  =  — . 

14 

Let  (h  :=  -I-  z£  denote  the  finite-element  approximation  in  the  uniform- 

patch  5(x°, Hza)  of  u  constructed  from  the  interpolant  and  the  periodic  finite- 
element  approximation  of  the  interpolation  error  (note  that  £h  can  be  constructed 
from  the  interpolant  Q™T  and  by  solving  a  small  discrete  problem  to  obtain  z£  in 
one  super-patch).  It  can  then  be  shown  that  ta  in  5(x°,  Hz0)- 

Theorem  2.  Let  the  assumptions  of  Theorem  1  hold  then 

lllt*fc  -  {‘  III*'**)  =  lll«*  -  QT  -  <  Cf*»HhHII^.^)  (61) 

where  C  is  independent  of  H  and  h 

In  (17)-  (44)  we  introduced  various  specific  element  error-indicators  i/r  which 
depend  on  the  finite  element  solution  uh  and  the  data-function  /.  The  error- 
estimate  in  Zo  is 
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£*=<*.(•**,/)  :=  {  £  (M«‘, /))’}’  • 

£z,  utilizes  the  finite  element  solution  uH  on  r  €  Th,  t  C  Z0,  and  on  the  neigh¬ 
borhood  of  Zo-  Let  us  denote  by  Z$  a  slightly  larger  subdomain  which  includes 
the  element  in  Z0  and  their  neighbors  (the  estimator  in  ZQ  depends  on  the  finite- 
element  solution  u*  in  Z$  and  the  function  /  (»«e  (9)))- 

We  will  impose  on  the  estimators  the  requirements  of  scale-invariance  and 
stability. 

Definition  1:  Let  X\  —  cx\ ,  =  cz 2  and 

z;  :={(*„  X,)  :  (*„*,)  €  Z.  }  ,  Ut(XuX,)  :=  ^(xux2)  , 

(62a) 

Tl  :={(X„.X,)  :  (*„*,)  €  Tk  }  ,  F')*,.*,)  :=  c~*f(x . 

Then  the  estimator  is  scale-invariant  if  for  any  c  6  (0,  oo) 

e*  (t5,f)-«*i(»‘./)  (6“) 

Remark  5.4.  The  energy-norm  is  invariant  with  respect  to  the  scaling  of  the  co¬ 
ordinates  and  hence  the  estimators  should  also  be  invariant.  On  the  other  hand  be¬ 
cause  the  differential  operator  is  of  second-order  the  scaling  appears  in  F(-Xi,-Xj) . 
It  is  obvious  that  every  error  estimator  should  be  scale-invariant. 

Definition  2:  The  estimator  Eg9  is  stable  if  for  any  Zo,  any  vh,  <ph  g  5^(Zq  )  and 
s,g  €  La(Z0)  and  0  <  a  <  1 

+  +  <c[a£j,(»‘,.)  +  (l  +  a-')fi,(/,J)]  (63) 

where  C  is  independent  of  h  and  the  class  of  admissible  meshes  under  consideration. 

Remark  5.5.  The  stability  of  the  error  estimator  describes  the  requirement  that 
small  changes  in  the  finite-element  solution  uh  and  the  data  /  should  result  in 
small  changes  in  the  value  of  the  estimator.  It  can  be  shown  that  all  estimators 
mentioned  above  are  scale-invariant  and  stable. 

The  following  theorem  states  that  whenever  an  estimator  is  stable  we  can  em¬ 
ploy  the  results  from  the  error-estimation  in  the  periodic  super-patch  to  study  the 
quality  of  the  estimator  in  the  actual  finite-element  mesh. 
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Theorem  3.  Let  the  assumptions  of  theorem  1  hold  and  the  estimator  be  scale- 
invariant  and  stable.  Farther  assume  that  for  the  estimators 


K(c(x»,*).^^)  =  ^M)>q>o  (M) 

(Here  we  use  the  notation  L'Q  :=  —A Q;  L'Q  is  a  constant  function.  Recall  that 
we  address  here  only  the  case  p  =  1  and  Q  is  a  polynomial  of  degree  2).  Then 

KI(S(x°,ffz.),<.\/)  =  K\S(x<‘,H2.),(k,L'Q)(l  +  Chi)  (65) 

Theorem  3  states  that  the  effectivity  index  K(S(nP,Hz0),uh,f)  is  asymptoti¬ 
cally  (as  h  — ►  0)  the  same  as  the  effectivity  index  #6(5(x?,JHj^),  £*,£'(?)  and  it  is 
easy  to  see  that 


*(S(x",ff*),«\£'G)  =  « (c(x°,h),('‘,L'Q)  (66) 

Hence  we  can  compute  the  effectivity  index  by  solving  a  small  discrete  problem 
in  one  super-patch  only  and  because  of  invariancy  with  respect  to  scalings  we  can 
study  only  the  master  mesh  on  the  master  super-patch.  Note  that  the  asympotical 
rate,  as  stated  in  the  theorem  3,  is  very  low.  Let  us  underline  that  the  aim  of 
the  theorem  was  to  show  the  asymptotic  behavior  and  not  to  obtain  the  optimal 
asymptotic  rate.  We  have  numerically  shown  (see  [53])  that  the  effectivity  index 
computed  from  the  cell  analysis  is  not  too  sensitive  to  the  surrounding  meshes  and 
that  the  term  Ch*  in  (65)  can  be  neglected. 

The  fact  that  we  have  used  square  super-patches  is  not  essential  (for  example 
it  is  possible  to  employ  regular  hexagons  or  parallelograms  as  super-patches).  For 
the  proofs  it  is  only  essential  that  the  mesh  is  locally  translation-invariant.  Let  us 
mention  that  if  the  solution  u  is  harmonic  (Au  =  0)  we  can  analyze  the  effectivity 
index  of  the  estimator  only  for  harmonic  polynomials.  Let  us  also  underline  that 
the  assumption  that  u  satisfies  Poisson’s  equation  and  the  elements  are  linear,  was 
made  only  for  simplicity.  The  theorems  hold  for  general  elliptic  operators  with  non¬ 
constant  coefficients  which  are  allowed  to  vary  smoothly  throughout  the  domain, 
and  for  higher-order  elements.  Although  the  theorems  were  formulated  for  the 
energy-norm  they  can  also  be  restated  for  other  norms  under  proper  assumptions. 


6  The  methodology  for  checking  the  estimators 

We  now  describe  the  numerical  methodology  which  is  employed  in  the  calculation 
of  the  robustness  index  for  a  given  estimator.  The  robustness  index  depends  on 
the  following  factors  (see  also  [50-53]): 
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a.  The  geometry  of  the  grid : 

The  main  factor  which  affects  any  error  estimator  is  the  geometry  of  the  grid 
namely  the  connectivity  or  topology  of  the  grid  and  the  geometry  ( distortion ) 
of  the  elements.  The  element  geometry  has  to  be  considered  in  connection  with 
the  differential  operator  for  the  boundary- value  problem.  For  example  the  or¬ 
thotropic  heat-conduction  operator  can  be  transformed  by  an  affine  transforma¬ 
tion  to  Laplace’s  operator  on  a  distorted  mesh  which  has  different  minimal  and 
maximal  angles. 

b.  The  class  of  solutions: 

In  this  paper  we  study  the  robustness  of  error  estimators  for  linear  elliptic 
equations  and  interior  mesh-subdomains.  It  is  well-known  that  the  solutions  of 
elliptic  boundary-value  problems  are  analytic  in  the  interior  of  the  domain  if  the 
coefficients  of  the  operator  and  the  right-hand  side  are  analytic;  for  this  reason  we 
will  consider  the  class  of  smooth  solutions.  Of  particular  interest  is  the  subclass 
of  smooth  functions  which  satisfies  the  homogeneous  differential  equation  (we  will 
refer  to  such  solutions  as  “harmonic”;  if  the  solution  satisfies  Laplace’s  equation  it 
is  truly  harmonic).  The  asymptotic  properties  of  error  estimators  for  the  class  of 
smooth  solutions  can  be  studied  by  considering  the  class  of  polynomials  of  degree 
(p+1). 

We  now  give  an  outline  of  the  numerical  procedure  which  determines  the  ro¬ 
bustness  index  TL  for  any  error  estimators. 

Let  us  assume  that  we  would  like  to  determine  the  robustness  of  an  estimator 
for  a  mesh  7*  from  a  given  class  of  meshes  (the  class  may  be  defined  as  the  class 
of  meshes  produced  by  a  commercial  mesh-generator  or  an  adaptive  code).  Let 
{u;-x}xUi  the  cells  of  elements  connected  to  the  vertices  of  the  mesh.  Let  a >x 
an  interior-cell  (it  is  separated  by  several  layers  of  elements  from  the  boundary, 
singular-points  and  material-interfaces).  Let 

u/o  :=  ufj 

and  for  s  >  1,  integer,  define  recursively 

u>*  :=  U  ux 

XeV(r') 

From  the  analysis  (given  in  [53]  for  locally-periodic  meshes)  and  numerical  experi¬ 
ence  (for  more  general  meshes)  we  know  that  the  effectivity  index  for  any  estimator 
in  the  cell  depends  practically  on  the  geometry  of  the  mesh  in  the  patch  u£  for 
only  a  small  s  (s  =  1  -  4)  (more  precisely  increasing  s  will  change  the  robustness 
index  only  minimally).  Given  a  mesh  7J,  (resp.  class  of  meshes  T)  constructed  by 
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a  mesh  generator,  there  are  patches  u>*  of  various  types  (topologies).  We  can  now 
analyze  all  these  patches  separately  and  compute  Tlrk  (resp.  H{T) )  as  defined  in 
(5b).  For  example,  in  Fig.  9  (resp.  Fig.  10)  we  indicate  the  mesh-cells  u>£  shaded 
and  the  patch  (resp.  u>J)  with  its  perigram  shown  by  thick  line,  for  various 
interior-cell/patch  combinations  from  the  grid  shown  in  Fig.  6  (resp.  Fig.  1).  We 
can  determine  the  robustness  of  any  error-estimator  in  Wg  as  follows: 

1.  Completion  of  the  mesh-cell  to  a  periodic  super-patch. 

We  translate  and  scale  the  subdomain  u/,  so  that  its  image  <l£  C  c  =  [—1, 1]  x 
[ — 1]  (see  Fig.  7b).  Then  we  employ  a  mesh-generator  to  complete  the  mesh  in 
into  a  periodic-mesh  T  in  the  super-patch  c  (see  Fig.  7c). 

2.  Periodic  boundary-value  problem: 

For  given  exact  solution  u  with  components  being  homogeneous  polynomials 
of  degree  (p  +  1)  denoted  by  Pp+i  (u  6  {Pp+ 1)2  or  exact  “harmonic”  solution 
u  €  {P p+i)J  such  that  i<(u)  =  0,  t  =  1,2),  material  properties  and  mesh-pattem 
(the  local  geometry  of  the  grid  in  the  cell  and  the  patch  )  determine 

a)  The  finite  element  solution  u/l , 

b)  The  exact  error  eh, 

by  solving  the  following  periodic  boundary-value  problem: 

Find  zh  6  (c)  such  that 

ftv.(*\vk)  =  Bn„,(u  -  u,\v»)  V  v*  6  H^„  ({!,„)  (67) 

where 

K-W  ={wk  6  H‘(2)  s  w*(-l,*t)  —  w*(l,na),  wh(a5lt— X)  =  wk(i,,  1) 
and  w»|r€( V,)\  ref}, 

u(*  denotes  a  continuous  piecewise  p-degree  interpolant  of  the  exact  solution  u. 

We  select  the  solution  zh  of  (66)  which  satisfies 

Jzh  =  J(u-  u?)  (68) 

i  i 

and  we  define 

uh  =  uf  +  z*  ;  e*  =  (u  -  uf )  -  zh  .  (69) 
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Thus  from  (69)  we  obtain  the  finite-element  solution  uh  and  the  exact  error  ef1  and 
for  a  given  estimator  we  compute  the  effectivity  index 


kuh  =  Kuh  (material  coeff.,  grid-material  orientation,  pattern,  solution  coeff.) 


Remark  6.1.  We  note  that  the  effectivity  index  k  is  independent  of  the  normal¬ 
ization  of  the  coefficients  in  (‘Pp+i)J. 

S.  Numerical  optimization 

We  let 


max  max  /c  >  , 

material  caaff.  solution  coeff.  0 


min 

material  caaff. 


miTi 

solution  caaff. 


(TO) 


The  bounds  C£° ,  C£°  for  a  given  class  of  solutions  and  materials,  a  given  pattern 
and  for  a  given  grid-material  orientation  can  now  be  computed  using  numerical 
optimization.  We  compute  the  robustness  of  the  estimator  in  wj  C  uj  from 


Huh  :=  max  Tl\  , 

"•  <=1,2,3,...  "O 


TV.h  :=  max 


{ (ii  -  of  i + u  -  cf  i) ,  (ii  -  4pi + n  -  -^0} 


(71) 


C?  C? 

We  can  then  compute  the  robustness  of  the  estimator  for  the  mesh  Th  by  computing 


'■=  max  TLuh  (72) 

“« —x  • 

Here  A/int(7h)  denotes  the  interior  vertices  of  the  mesh  (the  vertices  must  be  sep¬ 
arated  by  two  (or  more)  layers  of  elements  from  the  boundary  of  the  domain  and 
from  material-interfaces). 

Let  us  underline  the  reason  for  the  optimization.  In  general  we  know  only 
that  the  solution  satisfies  the  differential  equation,  e.g.  it  is  harmonic  (when  /  = 
0).  Hence  we  wish  that  the  effectivity  index  is  small  for  the  entire  considered 
class.  From  the  results  of  the  theoretical  study  (see  [53])  we  can  restrict  the  class 
of  functions  to  harmonic  polynomials  of  degree  (p  +  1)  only.  Establishing  the 
bounds  for  the  classes  of  meshes  and  solutions  of  interest  is  essential  for  judging 
the  robustness  of  an  estimator. 

Remark  6.2.  Note  that  the  methodology  assumed  that  the  local  mesh-size 
h  is  sufficiently  small  so  that  the  exact  solution  u  can  be  replaced  by  its  local 
Taylor  expansion  of  degree  (p+ 1)  (the  validity  of  the  methodology  is  asymptotic). 
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However  numerical  studies  show  that  the  robustness  index  governs  the  performance 
of  estimators  in  the  range  of  the  practical  engineering  computations  which  often 
employ  relatively  coarse  grids. 


7  Numerical  studies  of  robustness  of  various 
error-estimators 

We  present  examples  of  the  application  of  the  methodology  described  earlier  to 
study  the  robustness  of  several  error  estimators,  namely: 

1.  Implicit  element  residual  (Est.  1)  (Eqs.  (17)-(21))  [8],  [11],  [13],  [14]. 

2.  Implicit  element  residual  with  equilibration  (Est.  2)  (Eqs.  (17)-(36))  [8],  [17], 
[25],  [29],  [59],  [60],  [61]. 

3.  Implicit  element  residual  based  on  complementary  energy  principle  (Est.  3) 
(defined  below)  [25],  [26],  [27],  [28],  [32],  [60],  [61]. 

4.  Estimators  based  on  smoothening  or  Z-Z  estimators  (Est.  4)  (Eqs.  (37)-(44)) 
[40],  [41],  [42],  [43],  [44]. 

Note  that  in  [53]  we  also  studied  the  performance  of  several  other  estimators  (ex¬ 
plicit  element-residual,  subdomain  residual  etc.). 

We  now  proceed  to  the  discussion  of  the  numerical  studies. 

7.1  The  robustness  index  for  polynomial  solutions 

Here  we  address  three  questions: 

1)  How  many  layers  of  elements  should  be  taken  around  the  cell  a/£,  for  the 
values  of  Cya ,  C£° ,  71uk  to  be  practically  independent  of  the  surrounding 
mesh  in  the  periodic  super-patch  c  ? 

2)  How  much  do  C£°,  C%0  and  7iuh  computed  from  a  periodic  super-patch 
which  includes  w*  (s  =  3  for  meshes  of  triangles  and  s  =  4  for  meshes 
of  quadrilaterals)  differ  from  the  same  quantities  computed  from  the  entire 
original  mesh? 

3)  What  is  the  robustness  index  for  the  various  estimators? 


26 


We  will  consider  the  interior-patterns  (cell/patch  combinations)  shown  in  Figs.  10 
and  9  which  occur  in  the  meshes  shown  in  Figs.  1  and  6,  respectively. 

For  the  scalar  elliptic  problem  we  chose  the  exact  solution  to  be  either  a  general 
(homogeneous)  polynomial  of  degree  (p  +  1) 

QG{zux2)  =  Y^Oijxixi  ,  i  +  j=p+ 1,  t,j>0  (73) 


or  a  harmonic  polynomial  by  imposing  the  constraint 
example  for  Kut  =  &kt  we  have 


2 


Dk{KktDiQB)  =  0;  for 

M=i 


Qa(xi,z2)  =  ai(xl  -  x\)  +  a2ziz3,  forp+l=2  (74a) 


and 


QB{xux2)  =  al{x\-3xlx\)  +  a2{3x\x2-x\),  forp+l  =  3  (746) 

In  the  case  of  the  vector-valued  model  problem  (linear  elasticity)  we  employed 
homogeneous  “harmonic”  polynomial  solutions  of  degree  (p  +  1) 

Q  Be(fp+1)2:  ^D,K(  QH))=  0,  *=1,2  (75) 

3= * 

We  determined  the  robustness  index  of  the  error  estimators  for  the  cells 
using  the  approach  of  the  paper.  In  Table  1  we  give  the  values  of  <7£° ,  C(j° ,  7Lwh 
for  Est.  1  for  the  scalar  elliptic  problem  with  p  =  1,  for  the  pattern  1  (shown  in 
Fig.  9a),  when  s  =  1,  2,  3,  4,  5  layers  of  elements  around  are  taken  in  the 
patch  u>,  (see  Fig.  8).  In  Tables  2  and  3  we  give  the  values  of  C£a ,  C%° ,  and 
(s  =  3  for  the  meshes  of  triangles  and  s  =  4  for  the  meshes  of  quadrilaterals) 
for  the  Est.  1,  Est.  2,  and  Est.  4  for  the  scalar  elliptic  problem  with  p  =  1  and 
2  respectively.  These  are  compared  with  the  values  of  the  effect ivity  indices  C£° , 

Cy0  and  "R.  (%  is  obtained  by  substituting  ,  C$  into  (5a))  obtained  from  the 
finite-element  solution  of  a  Dirichlet  boundary-value  problem  in  the  grids  shown  in 
Figs.  1,  6  (the  big  mesh)  with  data  consistent  with  the  homogeneous  polynomial 

solutions  which  give  the  extremal  values  C£# ,  Cp°  of  the  effectivity  index  kuk  in 
the  periodic  super-patch  calculations.  The  periodic  super-patches  employed  for 
the  various  patterns  were  constructed  using  a  mesh-generator,  as  shown  in  Fig.  7. 
In  Table  4  a  similar  comparison  is  made  for  isotropic  elasticity  (Poisson’s  ratio 
v  =  0.3)  using  meshes  of  triangles. 
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From  the  reported  data  we  can  make  the  following  conclusions  related  to  the 
questions  formulated  above: 

1)  For  linear  elements  the  values  of  C^° ,  Cjj  ,  in  the  cell  u>£  do  not  practi¬ 
cally  change  when  the  patch  has  a  >  3  layers  of  elements  around  the  cell 
o/g .  For  higher  order  elements,  a  smaller  a  can  be  taken,  i.e.  a  =  1-2  (thus  for 
our  numerical  experiments,  we  considered  patches  with  a  =  3  for  the  meshes 
of  triangles  and  a  =  4  for  the  meshes  of  quadrilaterals). 

2)  For  a  given  exact  solution  which  is  a  homogeneous  polynomial  (harmonic  or 

uh  uh 

general)  of  degree  (p  +  1)  the  values  of  CL° ,  Cp  ,  obtained  from  a  periodic 

super-patch  are  very  close  to  the  values  Ci° ,  C^°  obtained  from  the  approx¬ 
imate  solution  of  a  boundary-value  problem  (with  data  compatible  with  the 
polynomial  solutions  which  correspond  to  the  extrema  of  the  effectivity  index 
in  the  periodic  super-patch)  obtained  using  the  big  mesh.  Note  that  we  did 
not  exactly  answer  question  2  because  we  did  not  perform  the  optimization 
in  the  big  mesh;  however  the  numerical  experience  from  [53]  indicates  that 

C£ ,  Cu°  are  very  close  to  the  extremal  values  of  the  effectivity  index  when 
the  optimization  is  performed  in  the  big  mesh. 

3a)  The  element-residual  without  equilibration  (Est.  1)  is  not  robust  for  the  gen¬ 
eral  meshes  (like  the  meshes  shown  in  Fig.  1  and  Fig.  6)  and  should  not  be 
used. 

3b)  The  equilibration  of  the  flux  increases  the  robustness  of  the  element-residual 
estimators. 

3c)  The  Z-Z  estimator  appears  to  be  the  most  robust. 

7.2  The  robustness  index  for  general  solutions 

In  7.1  above,  we  assumed  that  the  exact  solution  is  a  (homogeneous)  polynomial 
of  degree  (p  +  1) .  The  following  question  arises: 

Are  the  conclusions  based  on  the  assumption  that  the  solution  is  a  polynomial 
valid  for  general  solutions? 

We  now  give  an  example  which  shows  that  the  results  on  the  robustness  of 
estimators  obtained  using  polynomials  of  degree  (p  +  1)  give  a  good  indication  of 
the  local  performance  of  the  estimators  for  any  general  solution  (the  theoretical 
justification  for  patch  wise  uniform  meshes  is  given  in  theorem  3).  Let  us  consider 
the  general  solution 

u(*i,a:j)  =  (rx(zi,*j))*  sin(^0i(*i,zj))  +  (ra(*i,*a))*  sin(-0a(*i,*j)) 
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where 


**i(*»» *2)  =  [(*1  ~  +  (*a  “  3)s]  *  ,  *a)  =  tan-1  (^37)  » 

*  ®i  j 

r3(zuz2)  =  [(*!  —  3)a  -f  (*j  -  ^)J]a  ,  02(zuz2)  =  tan"1^*  — j)  , 

we  let  (2  =  (0,1)  x  (0,1)  and  let  7J,  be  the  grid  shown  in  Fig.  6.  We  solved 
the  Neumann  boundary-value  problem  in  Cl  using  (the  big  mesh)  and  data 
consistent  with  the  exact  solution  given  above  and  we  computed  the  effectivity 
index  for  Est.  1,  Est.  2  and  Est.  4  for  five  of  the  patterns  shown  in  Fig.  9.  We  also 
computed  the  effectivity  index  in  the  cells  by: 

a)  Solving  the  Neumann  problem  in  0  using  the  grid  7*  with  data  corresponding 
to  the  local  Taylor-expansion  (up  to  quadratic  degree)  about  the  central  node 
X  of  the  cell  =  u>x  • 

b)  By  using  the  local  Taylor-expansion  (about  the  central  node  of  the  patch) 
as  exact  solution  to  pose  the  periodic  boundary-value  problem  (66)  over  the 
periodic  super-patches  for  the  patterns  shown  in  Figs.  9a-9e  (patterns  1-5). 

The  effectivity  indices  for  the  cells  calculated  from  the  approximate  solution  (which 
was  obtained  using  the  three  different  ways  stated  above)  are  given  in  Table  5. 

We  observe  that  the  effectivity  indices  obtained  using  the  approximate  solu¬ 
tion  (computed  from  the  big  mesh  Th)  of  the  Neumann  boundary-value  problem 
formulated  using  data  obtained  either  from  the  general  solution  or  from  its  local 
Taylor  series  expansion  are  essentially  identical.  The  results  obtained  using  the 
local  Taylor  series  expansion  and  the  periodic  boundary-value  problem  (67)  are 
also  very  dose  to  the  values  obtained  from  the  other  two  types  of  problems;  still 
better  agreement  could  be  obtained  by  increasing  the  number  of  layers  «  in  wj 
(here  we  have  used  s  =  3  layers). 

7.3  The  influence  of  various  parameters  on  the  robustness 
of  the  estimators 

The  definitions  of  Est.  2,  Est.  3  and  Est.  4  include  the  choice  of  some  free 
parameters ,  namely: 

1)  The  choice  of  the  correction  of  the  edge  residuals  for  Est.  2  and  Est.  3. 
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2)  The  number  and  location  of  the  sampling  points  involved  in  the  definition  of 
Est.  4. 

3)  The  discrete  space  of  polynomials  used  in  the  patch-projections  or  the  so¬ 
lution  of  the  element-residual  problem  (we  may  use  the  bubble-space  or  the 
complete  space  of  polynomials  of  degree  (p+  1)  or  some  intermediate  space). 

Here  we  show  that  the  approach  can  be  used  for  the  study  of  the  dependence  of 
the  robustness  index  on  the  parameters  involved  in  the  definition  of  the  estimators 
or  to  find  the  most  robust  version  of  an  estimator  for  the  grids  of  interest.  We  also 
note  that  the  methodology  can  be  used  to  design  robust  versions  of  estimators  by 
considering  additional  parameters  (like  for  example  the  location  of  the  sampling 
points)  as  variables  and  by  performing  the  optimization  for  the  finite  number  of 
cdl/patch-pattems  which  occur  in  the  grids  of  interest. 

To  show  the  influence  of  the  various  parameters  we  will  consider,  in  the  numer¬ 
ical  studies  given  below,  a  set  of  distorted  meshes  in  the  periodic  super-patches 
shown  in  Figs.  11,  12.  (These  meshes  were  selected  to  include  non-standard  mesh- 
topologies  and  elements  with  large  maximum  angle  and  aspect  ratios.) 

A.  Z-Z  estimators 

The  robustness  of  Z-Z  estimators  depends  on  the  choice  of  the  sampling  points. 
The  following  choices  are  studied  for  the  quadrilaterals  with  p  =  1: 

(i)  One  sampling  point  at  the  center  of  the  element  ( Type  1 )  as  shown  in  Fig.  4d. 

(ii)  Four  sampling  points  along  the  axes  of  the  master-element  ( Type  2)  as  shown 
in  Fig.  4e.  These  sampling  points  correspond  to  the  one-dimensional  Gauss- 
Legendre  integration  points  of  order  two  along  the  coordinate-axes  of  the 
master-element . 

(iii)  Four  sampling  points  ( Type  3)  at  the  integration-points  of  the  (2  x  2)  Gauss- 
Legendre  product-rule  (shown  in  Fig.  4f). 

We  studied  the  performance  of  the  above  estimators  for  the  periodic-grids  shown 
in  Fig.  11;  here  the  cell  u>£  coincides  with  the  periodic  super-patch.  In  Table  6  we 

give  the  values  of  C%° ,  CyQ  and  TIuk  for  these  cells  when  the  solution  is  a  general 
polynomial  of  degree  2.  We  observe  that  the  Z-Z  estimator  of  Type  3  which  uses 
sampling  points  at  the  integration  points  of  the  2x2  Gauss  product-rule  shows 
the  best  robustness  for  these  meshes. 

For  the  triangles  we  studied  the  following  choices  of  sampling  points  shown  in 
Fig.  4a,  b,  c. 

(i)  One  sampling-point  located  at  the  centroid  of  the  triangle  (Fig.  4a);  this 
point  will  be  used  only  for  linear  elements. 


(ii)  Three  sampling-points  located  at  the  midsides  of  the  edges  of  the  triangle 
(Fig.  4b). 

(iii)  Seven  sampling  points  which  include  the  points  in  (i)  and  (ii)  and  the  mid¬ 
points  of  the  segments  which  connect  the  centroid  with  the  vertices  of  the 
triangle  (Fig.  4c). 

We  computed  the  upper-  and  lower-bound  of  the  effectivity-index  and  the  ro¬ 
bustness  index  for  the  periodic  meshes  shown  in  Figs.  12;  the  periodic  meshes  have 
given  aspect-ratio,  1/a  (see  Fig.  12a)  and  include  some  elements  with  very  large 
maximum  angle.  For  this  example  we  took  aspect-ratio  of  1/2  for  the  periodic 
meshes  shown  in  Figs.  12a-12d,  and  aspect-ratio  of  1/8  for  the  periodic  meshes 
shown  in  Figs.  12e,  12f.  In  Tables  7,  8  we  give  the  results  for  p  =  1  and  p  =  2, 
respectively,  obtained  by  solving  Laplace’s  equation.  In  Table  9  we  give  the  results 
for  p  =  1  for  isotropic  linear  elasticity.  It  should  be  noted  that  the  Z-Z  estimator 
is  exact  for  the  mesh  shown  in  Fig.  12e  (for  linear  elements).  We  observe  that  for 
the  meshes  shown  in  Fig.  12: 

(i)  For  Laplace’s  equation  and  p  =  1  the  choice  of  sampling  points  shown  in 
Fig.  4b  (3  sampling  points)  shows  the  best  robustness. 

(ii)  For  Laplace’s  equation  and  p  =  2  the  choice  of  sampling  points  shown  in 
Fig.  4c  is  the  most  robust. 

(iii)  For  isotropic  elasticity  and  p  =  1  the  choice  of  sampling  points  shown  in 
Fig.  4a  gives  the  most  robust  Z-Z  estimator. 

To  take  into  account  the  nature  of  the  exact  solution  (when  for  example  it 
is  harmonic)  we  may  consider  another  version  of  the  Z-Z  estimator  which  em¬ 
ploys  harmonic  polynomials  in  the  patch-projections.  To  check  if  harmonic  patch- 
projections  will  improve  the  robustness  of  the  Z-Z  estimator  for  grids  of  general 
quadrilaterals  we  employed  the  Z-Z  estimator  with  one-sampling  point  per  element 
and  the  periodic  cells  shown  in  Fig.  11.  The  results  for  the  values  of  C^° ,  Cffl , 
Huh  for  the  Z-Z  estimator  (for  the  Laplacian)  with  the  general  or  the  harmonic 
patch-projection  are  given  in  Table  10.  From  these  results  we  conclude  that  the 
Z-Z  estimator  is  less  robust  when  the  harmonic  patch-projection  is  employed. 

B.  Element  residual  estimators  v/iih  various  types  of  equilibration 

The  robustness  of  the  residual  estimators  depends  on  the  definition  of  the 
residual-functional.  Here  we  considered  two  types  of  equilibration  methods,  namely, 
Ladeveze’s  and  Ohtsubo-Kitamura’s  techniques.  Ladeveze’s  method  splits  the  in¬ 
terelement  fluxes  by  minimizing  the  functional  given  in  (31)  subject  to  the  con¬ 
straints  (29),  as  discussed  earlier.  We  considered  three  choices  of  the  weight  Wi  in 
the  definition  of  (31).  These  are: 
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(i)  Wi  =  1  (see  for  example  [25],  [32])  (equilibration  of  Type  A); 

(ii)  Wi  =  | c,- 1 “ 1  (see  [60],  [61])  (equilibration  of  Type  B); 

(iii)  Wi  =  |ej|_1,  where  c[  denote  the  length  of  side  e*  in  the  transformed,  plane  in 
which  the  operator  is  the  Laplaci&n  (equilibration  of  Type  C). 

We  computed  the  robustness  of  the  element-residual  estimator  with  Ladeveze’s 
equilibration  of  Type  A,  B,  C  for  the  periodic  meshes  shown  in  Figs.  11,  12a-12e. 
It  should  be  noted  that  the  estimator  is  exact  for  linear  elements  and  the  mesh 
shown  in  Fig.  12f,  while  for  quadratic  elements  it  is  exact  for  the  mesh  in  Fig.  12e. 
The  results  are  given  in  Tables  11-17.  Table  11  gives  the  results  for  bilinear  quadri¬ 
laterals  for  Est.  2  with  equilibrations  of  Type  A  and  C  and  for  the  periodic  meshes 
shown  in  Fig.  11.  Tables  12-15  give  the  results  for  the  meshes  of  triangular  elements 
shown  in  Fig.  12.  In  particular  Tables  12, 13  give  the  results  of  the  optimization  for 
Est.  2  ( Type  A  and  B)  for  linear  and  quadratic  elements,  respectively  (note  that 
we  selected  aspect-ratio  1/2  for  the  meshes  shown  in  Figs.  12a-12d  and  aspect- 
ratio  of  1/8  for  the  meshes  shown  in  Figs.  12e,  12f).  Table  14,  15  give  the  results 
of  the  optimization  for  the  model  problem  (9)  (orthotropic  heat-conduction)  with 
1  <  AT — / <  10  (the  ratio  K — IK-.~  was  included  as  a  parameter  in  the 
optimization)  and  various  orientations  of  the  principal  axes  of  orthotropy  with 
respect  to  the  periodic  meshes  shown  in  Fig.  12c,  Fig.  12e  respectively  (here  we 
selected  aspect-ratio  1/1).  Table  16  gives  the  results  for  the  robustness  of  Est.  2 
( Types  A  and  B)  for  the  model  problem  of  isotropic  elasticity  (with  Poisson’s  ratio 
v  =  0.3 ).  It  should  be  noted  that  for  the  elasticity  problem  we  used  aspect-ratio 
of  1/2  for  the  meshes  shown  in  Figs.  12a-  12d  and  aspect-ratio  of  1/8  for  the  mesh 
shown  in  Fig.  12e.  We  observe  that: 

(i)  For  the  meshes  of  quadrilaterals,  where  we  employed  Est.  2  (with  equilibra¬ 
tion  of  Type  A  and  C),  the  robustness  depends  on  the  mesh-topology;  both 
types  have  reasonable  robustness. 

(ii)  For  the  meshes  of  triangles  Est.  2  with  equilibration  of  Type  B  has  the  best 
robustness  for  the  Laplacian. 

(iii)  From  the  results  for  the  case  of  orthotropic  heat-conduction  we  cannot  make 
a  general  statement  about  which  type  of  equilibration  is  the  best.  All  three 
types  seem  to  give  reasonably  robust  element-residual  estimator. 

(iv)  For  the  case  of  isotropic  elasticity  the  equilibration  of  Type  B  produces  the 
most  robust  estimator. 

We  also  implemented  element-residual  estimators  with  Ohtsubo-Kitamura  flux- 
splitting  for  the  isotropic  elasticity  problem.  The  element-residual  problem  is  de¬ 
fined  as  in  (21)  with  the  right-hand  side  defined  (with  or  without  equilibration)  as 
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described  in  3.2.2.  The  discrete  space  employed  in  the  solution  of  the  local  problem 
was  defined  as: 

(i)  The  discrete  space  which  corresponds  to  the  8-node  serendipity  element  ( com¬ 
plete  serendipity  space). 

(ii)  The  discrete  space  spanned  by  the  shape-functions  of  the  mid-side  nodes  in 
the  8-node  serendipity  element  ( bubble  serendipity  space). 

In  Table  17  we  give  the  values  of  C£° ,  Cu° ,  for  the  element-residual  estima¬ 
tors  with  Ohtsubo-Kitamura’s  flux-splitting  (isotropic  elasticity  with  v  =  0.3).  In 
columns  2-4  we  give  the  results  when  th*  ;tial  flux-splitting  (eq.  (33))  is  employed 
(without  equilibration)  together  with  the  rendipity  bubble-space.  In  columns  5-7 
(resp.  8-10)  we  give  the  results  when  the  equilibration  is  employed  together  with 
the  complete  (resp.  bubble)  serendipity  space.  We  observe  the  following: 

(i)  When  the  complete  serendipity  space  is  employed  (together  with  equilibra¬ 
tion)  the  estimator  lacks  robustness. 

(ii)  When  the  bubble  serendipity  space  is  used  together  with  equilibration  a 
robust  estimator  is  obtained. 

C.  Estimators  based  on  the  complementary  energy  principle 

Error  estimators  based  on  the  complementary  energy  principle  (e.g.  [25],  [26], 
[27],  [28],  [60],  [61])  are  often  preferred  by  engineers  because  they  provide  guaran¬ 
teed  upper-bounds  of  the  error.  Here  we  use  the  methodology  of  the  paper  to  check 
the  robustness  of  a-posteriori  error  estimators  based  on  the  complementary-energy 
principle  for  the  orthotropic  heat-conduction  and  elasticity  problems  and  the  class 
of  periodic  grids  shown  in  Figs.  11, 12.  The  error  estimators  are  defined  as  follows: 

a.  Orthotropic  heat-conduction: 


Vr  := 


where 


(76) 


i  B‘r(e =  *«(»,)  V.,E TVf,  (77) 


b.  Elasticity 


ifr  :=  ||6-ff(u»)||i»w,.-i 


(78) 


where  a  €  JMj  such  that 
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-V-<T  =  f 


in  r 


V  c  6  E{r) 


(79) 


Here  JMj  denotes  the  stress-space  for  the  Johnson- Merrier  composite  equilibrium 
triangular  element  (see  [63]).  (We  employed  the  explicit  solution  of  (79)  given  in 
[601) 

For  the  problem  of  orthotropic  heat-conduction  we  computed  the  robustness 
of  Est.  3  using  Ladeveze’s  flux-splitting  with  Type  A  and  Type  B  equilibrations; 
the  results  are  given  in  Table  18.  For  the  elasticity  problem  we  repeated  the 
computations  for  the  above  types  of  equilibration  and  for  the  periodic  meshes 
shown  in  Fig.  12a,  b,  c;  the  results  are  given  in  Table  19.  It  should  be  noted  that 
we  employed  aspect-ratio  of  1/2  for  the  meshes  in  Fig.  12a-12d  and  aspect-ratio  of 
1/8  for  the  mesh  shown  in  Fig.  12e.  We  observe  that: 


a)  For  both  cases  (scalar  and  vector-valued  problem)  the  estimator  with  equili¬ 
bration  of  Type  B  shows  superior  robustness  than  the  estimator  with  equili¬ 
bration  of  Type  A. 

b)  The  error  estimators  based  on  the  complementary  energy  principle  are  less 
robust  (by  far)  than  the  corresponding  element-residual  estimators. 


8  Summary  of  conclusions 

1.  The  validation  of  the  performance  of  the  estimators  based  on  the  robustness 
index  allows  objective  comparisons  between  the  various  error  estimators. 

2.  A  numerical  methodology  for  computing  the  robustness  index  is  given. 

3.  The  methodology  takes  directly  into  account  the  factors  which  affect  the 
performance  of  estimators  namely  the  geometry  of  the  grid,  the  differential 
operator  and  the  nature  of  the  solution.  The  methodology  has  theoretical 
basis  and  can  be  used  to  study  the  robustness  of  error  estimators  for  the 
complex  grids  which  are  used  in  engineering  computations. 

4.  The  methodology  allows  us  to  check  the  quality  of  any  new  estimator  even  if 
it  is  only  available  as  a  black-box  computer  subroutine. 

5.  It  is  possible  to  use  the  methodology  to  maximize  the  robustness  of  a  given 
estimator  for  a  class  of  meshes  of  interest. 
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6.  The  methodology  addresses  the  robustness  index  for  the  estimators  for  el¬ 
ements  in  the  interior  of  the  domain  and  smooth  solutions.  The  case  of 
unsmooth  solutions  and  elements  at  the  boundary  will  be  addressed  in  a 
forthcoming  paper. 

7.  The  Z-Z  estimator  seems  to  be  the  most  robust.  For  the  bilinear  quadrilat¬ 
eral  elements  the  sampling  points  defined  in  Fig.  4f  are  recommended.  For 
the  linear  (resp.  quadratic)  triangular  elements  the  choice  shown  in  Fig.  4b 
(resp.  4a)  is  the  best  for  Laplace’s  equation,  while  for  linear  triangles  and 
elasticity  the  choice  4a  is  the  most  robust.  Thus,  for  linear  elements  the 
choice  4a  or  4b  could  be  recommended. 

8.  The  element  residual  estimators  should  be  used  only  with  equilibration. 
Ladeveze’s  equilibration  of  Type  B  is  recommended. 

9.  The  Ohtsubo-Kitamura  flux-splitting  technique  gives  a  robust  estimator  only 
when  a  bubble  space  is  employed  in  the  element-residual  problem. 

10.  The  estimators  based  on  the  complementary  energy-principle  have  poor  ro¬ 
bustness  and  are  not  recommended. 

We  remark  that  the  conclusions  made  above  are  related  to  the  use  of  general 
meshes.  Use  of  families  of  particular  meshes  could  influence  the  conclusions. 
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Table  1.  Influence  of  the  size  of  the  patch  on  the  value  of  the  robustness- 
index  obtained  from  the  periodic  super-patch.  Laplace’s  equation,  quadratic 
harmonic  polynomial  solution,  linear  elements  (p  —  1).  Pattern  1  (shown  in  Fig.  9a) 
with  s  =  1,  2,  3,  4,  5  layers  around  the  cell  (as  shown  in  Fig.  8)  is  employed  in 
the  computation  of  the  robustness  for  Est.  1. 
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Table  2a.  Accuracy  of  the  methodology  for  general  meshes:  Laplace’s  equation, 
quadratic  harmonic  polynomial  solution,  linear  elements  (jp  —  1).  The  mesh  of 
triangles  shown  in  Fig.  6  and  the  patterns  1-7  shown  in  Figs.  9a-9g  are  employed 
in  the  computation  of  the  robustness. 
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Table  2b.  Accuracy  of  the  methodology  for  general  meshes:  Laplace’s  equation, 
quadratic  harmonic  polynomial  solution,  linear  dements  (p  =  1).  The  mesh  of 
quadrilaterals  shown  in  Fig.  1  and  the  patterns  1-5  shown  in  Figs.  lOa-lOe  are 
employed  in  the  computation  of  the  robustness. 
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Table  3.  Accuracy  of  the  methodology  for  general  meshes:  Laplace’s  equation, 
cubic  harmonic  polynomial  solution,  quadratic  elements  (p  =  2).  (a)  The  mesh  of 
triangles  shown  in  Fig.  6  and  the  patterns  1-7  shown  in  Figs.  9a-9g  are  employed 
in  the  computation  of  the  robustness. 
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Isotropic  Elasticity,  Linear  Triangles 

Element  residual  with  equilibration  (Eat.  2) 
(Ladeveze’g  equilibration,  eq.  (32)) 


ZZ  estimator  (Eat.  4) 


Table  4.  Accuracy  of  the  methodology  for  general  meshes:  Isotropic  elasticity, 
quadratic  “harmonic”  polynomial  solution,  linear  elements  (p  =  1).  The  mesh  of 
triangles  shown  in  Fig.  6  and  the  patterns  1,  3,  4,  5  shown  in  Fig.  9a,  c,  d,  e  are 
employed  in  the  computation  of  the  robustness. 


Pattern 


Periodic  Problem 


1.054  0.062 


1.002 

1.043 


Dirichlet  BVP 


1.027  |  0.035 
0.051 
0.997  0.077 
.946  I  1.002  0.059 


Dirichlet  BVP 


Of  0$  K 


.899  1.005 
.856  1.190 


Periodic  Problem 

Pattern 

cf 

cf 

1 

.965 

1.033 

0.068 

3 

g? 

.959 

1.049 

OiKil 

4 

Bj 

.875 

1.037 

0.179 

5 

1 

.853 

1.195 

0.342 
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Laplace’s  Equation,  linear  Triangles 


Pattern 

Neumann  BVP 

Periodic  Problem 

General  solution 

Taylor  series 

Taylor  series 

K 

(Est.  1) 

K 

(Est.  2) 

K 

(Est.  3) 

K 

(Est.  1) 

K 

(Est.  2) 

K 

(Est.  3) 

K 

(Est.  1) 

1C 

(Est.  2) 

K 

(Est.  3) 

1 

1.219 

1.220 

0.993 

0.999 

1.013 

2 

1.293 

1 

1.293 

■ 

;  '  : 

1.298 

m wm 

1.009 

3 

1.007 

1.015 

1.006 

1.014 

0.991 

4 

1.003 

i 

1.086 

1.003 

1.098 

1.015 

1.003 

5 

m 

1.031 

0.983 

1.317 

1.031 

0.983 

1.316 

1.034 

0.987 

Table  5.  Applicability  of  the  methodology  for  general  solutions:  Laplace’s  equa¬ 
tion,  harmonic  solution,  linear  elements  (p  =  1).  Comparison  of  the  values  of 
the  effectivity  index  for  the  cells  computed  using  three  different  approximate  so¬ 
lutions:  The  solution  of  a  Neumann  boundary-value  problem  in  the  domain  and 
mesh  shown  in  Fig.  6  with  the  data  taken  from  the  exact  solution  (Columns  2,  3, 
4);  the  solution  of  a  Neumann  boundary- value  problem  in  the  domain  and  mesh 
shown  in  Fig.  6  with  the  data  taken  from  the  quadratic  Taylor-series  expansion 
about  the  central-node  of  the  mesh-cell  Ug  for  which  the  effectivity  index  is  com¬ 
puted  (Columns  5,  6,  7);  the  solution  of  a  periodic  boundary-value  problem  in 
super-patches  which  include  the  patches  shown  in  Figs.  9a-9e  with  the  data  taken 
from  the  quadratic  Taylor-series  expansion  about  the  central  node  of  each  of  the 
mesh-cells  (Columns  8,  9, 10). 


46 


Poisson’s  Equation,  Bilinear  Quadrilaterals 
Z-Z  Estimators 
Periodic  Mesh  1 


Table  6a.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z- 
Z  estimators,  Poisson’s  equation,  general  quadratic  polynomial  solution,  bilinear 
quadrilaterals  (p  =  1).  Comparison  of  the  robustness  of  the  various  versions  of  the 
Z-Z- estimator.  The  robustness  results  are  given  for  the  periodic  mesh  shown  in 
Fig.  11a. 


Poisson’s  Equation,  Bilinear  Quadrilaterals 

Z-Z  Estimators 

Periodic  Mesh  2 

Type 

Kvnm-r 

Kmm 

CL 

Cv 

n 

1 

1 

0.945 

0.073 

2 

1 

3 

1 

0.995 

1.009 

0.014 

1 

100 

0.860 

1.343 

0.483 

2 

100 

0.891 

0.438 

3 

100 

0.946 

0.239 

Table  6b.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z- 
Z  estimators,  Poisson’s  equation,  general  quadratic  polynomial  solution,  bilinear 
quadrilaterals  (p  =  1).  Comparison  of  the  robustness  of  the  various  versions  of  the 
Z-Z-estimator.  The  robustness  results  are  given  for  the  periodic  mesh  shown  in 
Fig.  lib. 


Poisson’s  Equation,  Bilinear  Quadrilaterals 

Z-Z  Estimators 

Periodic  Mesh  3 

Type 

III 

CL 

Cv 

n 

1 

1 

0.995 

1 

2 

1 

1.001 

3 

1 

0.999 

1 

mm 

1 

100 

0.956 

1.015 

2 

100 

0.955 

1.011 

3 

100 

0.989 

1.023 

0.034 

Table  6c.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z- 
Z  estimators,  Poisson’s  equation,  general  quadratic  polynomial  solution,  bilinear 
quadrilaterals  (p  =  1).  Comparison  of  the  robustness  of  the  various  versions  of  the 
Z-Z-estimator.  The  robustness  results  are  given  for  the  periodic  mesh  shown  in 
Fig.  11c. 
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Poisson’s  Equation,  Bilinear  Quadrilaterals 

Z-Z  Estimators 

Periodic  Mesh  4 

Type 

AjqIX 

Amin 

CL 

Cv 

n 

1 

l 

0.991 

2 

l 

m 

3 

l 

0.994 

1.005 

1 

100 

0.849 

1.021 

p  i 

2 

100 

0.848 

3 

100 

0.847 

1 

0.190 

Table  fid.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z- 
Z  estimators,  Poisson’s  equation,  general  quadratic  polynomial  solution,  bilinear 
quadrilaterals  (p  =  1).  Comparison  of  the  robustness  of  the  various  versions  of  the 
Z-Z-estimator.  The  robustness  results  are  given  for  the  periodic  mesh  shown  in 
Fig.  lid. 
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Laplace’s  Equation,  Linear  Triangles 


Table  7.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z- 
Z-estimators,  Laplace’s  equation,  harmonic  quadratic  polynomial  solution,  linear 
triangles.  Comparison  of  the  robustness  of  the  various  versions  of  the  Z-Z  estima¬ 
tor.  The  robustness  results  are  given  for  the  periodic  meshes  shown  in  Fig.  12. 


Laplace’s  Equation,  Quadratic  Triangles 

Periodic 

ZZ  estimator  ( Type  2) 

ZZ  estimator  ( Type  S) 

_ 

Mesh 

CL 

Cu 

n 

CL 

Cu 

n 

1 

m 

0.049 

■ 

BH 

0.027 

2 

0.893 

EES! 

0.122 

Ha 

1 

3 

0.782 

0.990 

W:I'] 

0.873 

0.174 

4 

0.823 

1.003 

0.843 

'  *■ 

0.192 

6 

1.001 

1.006 

0.007 

1.013 

1.047 

0.060 

Table  8.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z- 
Z-estimators,  Laplace’s  equ&.ion,  harmonic  cubic  polynomial  solution,  quadratic 
triangles  (p  =  2).  Comparison  of  the  robustness  of  the  various  versions  of  the 
Z-Z  estimator.  The  robustness  results  are  given  for  the  periodic  meshes  shown  in 
Fig.  12. 
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Isotropic  Elasticity,  Linear  Triangles 


Periodic  ^  estimator  (Type  1)  Z Z  estimator  (Type  S)  ZZ  estimator  (Type  S) 

Mesh  CL  Cv  H  CL  \  Cv  \  K  CL  Cv  H 


0.944 
0.969 
0.726  1.044 

0.693  1.039 

00 


Table  0.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z-Z- 
estimators,  isotropic  elasticity  problem,  “harmonic”  quadratic  polynomial  solution, 
linear  triangles  (p  =  1).  Comparison  of  the  robustness  of  the  various  versions  of 
the  Z-Z  estimator.  The  robustness  results  are  given  for  the  periodic  meshes  shown 
in  Fig.  12. 
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Laplace’s  Equation,  Bilinear  Quadrilaterals 

Z-Z  Estimators 

Periodic 

Harmonic  projection 

Original  projection 

Mesh 

CL 

Cu 

K 

Cl 

Cu 

K 

1 

■ 

1.442 

0.444 

0.999 

B 

HEl 

2 

1.551 

0.553 

tm 

3 

1.465 

0.465 

1.017 

0.022 

4 

0.999 

1.301 

0.302 

0.999 

1.012 

■Hi 

Table  10.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Z-Z 
estimators,  Laplace’s  equation,  bilinear  quadrilaterals.  Comparison  of  the  robust¬ 
ness  of  two  versions  of  the  Z-Z  estimator  based  on  the  general  or  the  harmonic 
patch  projection  for  the  periodic  meshes  shown  in  Fig.  11. 
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Table  11.  Influence  of  various  parameters  on  the  robustness  of  estimators: 
Element-residual  error-estimators  with  various  types  of  equilibration.  Laplace’s 
equation,  quadratic  harmonic  polynomial  solution,  bilinear  quadrilaterals  (p  =  1). 
Comparison  of  the  robustness  of  estimators  Est.  2/  Type  A ,  Est.  2/  Type  C  for  the 
periodic  meshes  shown  in  Fig.  11. 
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Laplace’s  Equation,  Linear  Triangles 


Element  residual  with  equilibration 

Periodic 

Type  A  (in,- 

=  1) 

Type  B  (tu<  = 

W1) 

Mesh 

CL 

Ci r 

n 

CL 

Cv 

n 

1 

0.958 

1.002 

0.957 

PH 

2 

1.226 

j 

:  ;  -  ^ 

BBjHi 

3 

1.636 

4 

HRS 

1.266 

Ht%fi 

0.113 

5 

H 

2.958 

1.958 

Q 

1.003 

0.003 

Table  12.  Influence  of  various  parameters  on  the  robustness  of  estimators: 
Element-residual  error-estimators  with  various  types  of  equilibration.  Laplace’s 
equation,  quadratic  harmonic  polynomial  solution,  linear  triangles  (p  =  1).  Com¬ 
parison  of  the  robustness  of  estimators  Est.  2/ Type  A,  Est.  2 /Type  B  for  the 
periodic  meshes  shown  in  Fig.  12. 
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Laplace’s  Equation,  Quadratic  Triangles 


Element  residual  with  equilibration 


Periodic 

Type  A  (wi 

=  1) 

Type  B  (to,-  = 

M-1) 

Mesh 

CL 

Cv 

K 

CL 

Cv 

n 

1 

0.944 

1.016 

0.075 

0.958 

0.984 

0.060 

2 

0.991 

1.020 

0.029 

0.991 

1.007 

0.016 

3 

1.007 

1.262 

0.269 

0.989 

1.003 

0.014 

4 

0.999 

1.082 

0.083 

0.919 

0.994 

0.094 

6 

1.000 

2.369 

1.369 

1.000 

1.001 

0.001 

Table  IS.  Influence  of  various  parameters  on  the  robustness  of  estimators: 
Element-residual  error-estimators  with  various  types  of  equilibration.  Laplace’s 
equation,  cubic  harmonic  polynomial  solution,  quadratic  triangles  (p  =  2).  Com¬ 
parison  of  the  robustness  of  estimators  Est.  2/ Type  A ,  Est.  2/  Type  B  for  the 
periodic  meshes  shown  in  Fig.  12. 
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Orthotropic  Heat-Conduction,  Linear  Triangles 


Element  residual  with  equilibration 


Grid-Material  Type  A  {wi  =  1)  Type  B  (ti =  je^-1!)  Type  C  ( vii  =  |ei|-1) 
Orientation 


Cu 

% 

CL 

Cu 

2.280 

0.962 

1.474 

1.727 

0.750 

1.002 

1.495 

1.638 

0.679 

0.995 

1.992 

1.693 

0.699 

0.994 

2.245 

1.435 

0.463 

0.994 

1.451 

1.948 

0.976 

0.960 

1.061 

2.280 

1.284 

0.962 

1.474 

1.727 

0.750 

0.995 

1.495 

1.638 

0.679 

0.995 

1.992 

1.693 

0.699 

0.994 

2.245 

1.435 

0.463 

0.994 

1.451 

1.948 

0.976 

0.960 

1.061 

2.280 

1.284 

0.962 

1.474 

.457 
.101  0.894 

.512  0.877 

.500  0.941 
0.995 
0.995 
0.993 
0.894 
0.877 


Cu 

K 

1.237 

m 

1.301 

0.360 

1.989 

0.994 

2.124 

1.129 

1.331 

0.338 

1.079 

0.192 

1.237 

0.360 

1.301 

0.360 

1.989 

0.994 

2.124 

1.129 

1.331 

0.338 

1.079 

0.192 

1.237 

0.360 

Table  14.  Influence  of  various  parameters  on  the  robustness  of  estimators: 
Element-residual  error-estimators  with  various  types  of  equilibration.  Orthotropic 
heat-conduction  problem  (1  <  K — /  Km m  <  10),  quadratic  “harmonic”  polyno¬ 
mial  solution,  linear  triangles  (p  =  1).  Comparison  of  the  robustness  of  estimators 
Est.  2 1  Type  A ,  Est.  2/  Type  B,  Est.  2/  Type  C  for  the  periodic  mesh  shown  in 
Fig.  12c  and  various  grid-material  orientations. 
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Orthotropic  Heat-Conduction,  Linear  Triangles 


Element  residual  with  equilibration 

Grid-Material 

Type  A  (in,-  = 

=  1) 

Type  C  = 

141 -1) 

Orientation 

* (deg) 

CL 

Cv 

n 

CL 

Cv 

n 

-90 

1.000 

1.419 

0.419 

1.000 

1.010 

0.010 

-75 

1.000 

1.156 

0.156 

0.968 

1.002 

0.035 

-60 

0.894 

1.000 

0.119 

0.871 

1.013 

0.161 

-45 

0.816 

1.000 

0.225 

0.816 

1.030 

0.255 

-30 

0.894 

1.000 

0.119 

0.871 

1.000 

0.148 

-15 

1.000 

1.156 

0.156 

0.968 

1.002 

0.035 

00 

1.000 

1.419 

0.419 

1.000 

1.010 

0.010 

15 

1.000 

1.156 

0.156 

0.968 

1.002 

0.035 

30 

0.894 

1.000 

0.119 

0.871 

1.000 

0.148 

45 

0.816 

1.000 

0.225 

0.816 

1.030 

0.255 

60 

0.894 

1.000 

0.119 

0.871 

1.013 

0.161 

75 

1.000 

1.156 

0.156 

0.968 

1.002 

0.035 

90 

1.000 

1.419 

0.419 

1.000 

1.010 

0.010 

Table  15.  Influence  of  various  parameters  on  the  robustness  of  estimators: 
Element-residual  error-estimators  with  various  types  of  equilibration.  Orthotropic 
heat-conduction  problem  (1  <  K — /  <  10),  quadratic  “harmonic”  polyno¬ 

mial  solution,  linear  triangles  (p  =  1).  Comparison  of  the  robustness  of  estimators 
Est.  2/  Type  A ,  Est.  2/  Type  C  for  the  periodic  mesh  shown  in  Fig.  12f  and  various 
grid-material  orientations. 


Isotropic  Elasticity,  Linear  Triangles 

Element  residual  with  equilibration 

Periodic 

Mesh 

Type  A  ( Wi  =  1) 

Ttpt  B  (w,  =  |<|->) 

CL 

Cv 

n 

CL 

Ci j 

n 

1 

2 

3 

4 

5 

0.909 

0.991 

1.053 

0.965 

1.000 

1.052 

1.526 

2.526 
1.510 
5.292 

0.150 

0.535 

1.579 

0.545 

4.292 

0.908 

0.999 

1.007 

0.910 

0.996 

1.075 

1.053 

1.365 

1.374 

1.028 

0.171 

0.054 

0.372 

0.464 

0.032 

Table  16.  Influence  of  various  parameters  on  the  robustness  of  estimators: 
Element-residual  error-estimators  with  various  types  of  equilibration.  Isotropic 
elasticity,  quadratic  “harmonic”  polynomial  solution,  linear  triangles  (p  =  1). 
Comparison  of  the  robustness  of  estimators  Est.  2/ Type  A,  Est.  2/ Type  B  for 
the  periodic  meshes  shown  in  Fig.  12. 
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Table  17.  Influence  of  various  parameters  on  the  robustness  of  estimators: 
Isotropic  elasticity,  quadratic  “harmonic”  polynomial  solution,  bilinear  quadrilat¬ 
erals  (p  =  1).  Element-residual  error-estimator  without  equilibration  using  the 
residual  functional  in  (33)  (columns  2-4)  and  with  Ohtsubo-Kitamura  equilibra¬ 
tion  (using  the  complete  serendipity  space  (columns  5-7)  and  using  the  serendipity 
bubble  space  (columns  8-10)).  Comparison  of  the  robustness  of  the  estimators  for 
the  periodic  meshes  shown  in  Figs,  lla-llf. 
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Laplace’s  Equation,  Linear  Triangles 

Estimators  based 

on  complementary  energy 

Periodic 

Equilibration  of  Type  A 

Equilibration  of  Type  B 

Mesh 

CL 

Cv 

n 

CL 

Cv 

n 

1 

1.520 

1.636 

1.156 

1.636 

1.156 

2 

1.254 

1.966 

1.220 

1.772 

1.026 

3 

1.641 

2.673 

2.314 

1.362 

1.742 

1.104 

4 

2.007 

4.205 

4.212 

1.792 

3.696 

3.488 

5 

1.225 

3.825 

3.050 

1.225 

1.735 

0.960 

Table  18.  Influence  of  various  parameters  on  the  robustness  of  estimator? -  Estima¬ 
tors  based  on  the  complementary  energy-principle.  Laplace’s  equation,  quadratic 
harmonic  polynomial  solution,  linear  triangles  (p  =  1).  Comparison  of  the  robust¬ 
ness  of  estimators  Est.  3/  Type  A  and  Est.  3/  Type  B  for  the  periodic  meshes  shown 
in  Fig.  12. 
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Isotropic  Elasticity,  Linear  Triangles 


Estimators  based 

on  complementary  energy 

Periodic 

Mesh 

Equilibration  of  Type  A 

Equilibration  of  Type  B 

CL 

Cl T 

% 

CL 

Cv 

n 

1 

1.184 

4.113 

3.297 

4.046 

2 

1.399 

7.579 

1 

5.902 

3 

14.399 

89.634 

102.033 

1 

65.359 

Table  10.  Influence  of  various  parameters  on  the  robustness  of  estimators:  Esti¬ 
mators  based  on  the  complementary  energy- principle.  Isotropic  elasticity  problem, 
quadratic  “harmonic”  polynomial  solution,  linear  triangles  (p  =  1).  Comparison 
of  the  robustness  of  estimators  Est.  3/  Type  A  and  Est.  3/  Type  B  for  the  periodic 
meshes  shown  in  Fig.  12a,  b,  c. 
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List  of  Figures 

Figure  1.  The  mesh-cell  u>£  (dark  gray)  and  the  surrounding  layers  of  elements 
(light  gray)  which  influence  the  error  (and  the  error-estimator)  in  u>£. 

Figure  2.  Edge  t  with  its  normal  u  and  the  elements  rw,  rta  connected  to  it. 

Figure  3.  (a)  The  vertex  X  with  the  elements  r*  and  the  edges  attached  to  it; 
(b)  The  local  enumeration  for  the  degrees  of  freedom  of  the  correction  6  for  edges 
connected  to  node  X. 

Figure  4.  Various  selections  of  sampling  points  in  the  master-triangle  and  rect¬ 
angle. 

Figure  5.  (a)  A  domain  fl  with  the  subdomains  S(x°,  i/°),  Q0,  Z0  which  are 
exactly  covered  by  a  periodic  array  of  super-patches;  (b)  An  example  of  finite- 
element  grid  made  by  the  periodic  repetition  of  a  periodic  super-patch;  (c)  The 
periodic  super-patch. 

Figure  6.  Typical  example  of  a  general  finite-element  grid  of  triangles  generated 
by  a  commercial  mesh-generator. 

Figure  7.  Extraction  of  a  patch  and  completion  to  a  periodic  super-patch  for 
meshes  of  triangles  and  quadrilaterals  (a)  The  actual  grid  of  triangles  with  the 
subdomains  u>£;  (b)  The  subdomain  w,  with  u;£  in  its  interior;  (c)  The  subdo¬ 
main  u>3  embedded  into  a  periodic  super-patch;  (d)  The  actual  grid  of  quadrilaterals 
with  the  subdomains  wj;  (e)  The  subdomain  u/*  with  u/£  in  its  interior;  (f)  The 
subdomain  w*  embedded  into  a  periodic  super-patch. 

Figure  8.  Influence  of  the  size  of  the  patch  on  the  calculation  of  the  robustness- 
index  7luk:  The  cell  ur£  (shown  without  shading)  surrounded  by  several  mesh-layers 
(indicated  by  various  tones  of  gray-shading). 

Figure  0.  General  mesh  of  triangles  generated  by  a  commercial  mesh-generator 
(shown  in  Fig.  6):  (a)-(g)  Cell/patch  combinations  (patterns)  1-7.  The  cell  u>£  is 
shaded  gray;  the  perigram  of  the  patch  is  shown  in  thick  black  line. 

Figure  10.  General  unstructured  mesh  of  quadrilaterals  (shown  in  Fig.  1):  Cell/ 
patch  combinations  (patterns)  used  in  the  study  of  the  robustness  index  for  the 
various  estimators.  The  cell  consists  of:  (a)  3  elements,  (b)  4  elements,  (c)  6 
elements,  (d)  7  elements,  (e)  8  elements  connected  to  a  node.  The  cell  tvj  is  shaded 
gray;  the  perigram  of  the  patch  u>*  is  shown  in  thick  black  line. 

Figure  11.  Influence  of  various  parameters  on  the  robustness  of  the  estimators: 
Periodic  meshes  of  bilinear  quadrilateral  elements,  (a)  Mesh  1;  (b)  Mesh  2;  (c) 
Mesh  3;  (d)  Mesh  4. 

Figure  12.  Influence  of  various  parameters  on  the  robustness  of  the  estimators: 
Periodic  meshes  of  triangles,  (a)  Mesh  1;  (b)  Mesh  2;  (c)  Mesh  3;  (d)  Mesh  4;  (e) 
Mesh  5;  (f)  Mesh  6. 
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